| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module netcdf_results_output | ||
| 2 | !> NetCDF output for particle tracing results. | ||
| 3 | !> | ||
| 4 | !> Provides structured binary output (results.nc) containing: | ||
| 5 | !> - Per-particle loss data (times_lost, trap_par, perp_inv, classification) | ||
| 6 | !> - Phase space positions (zstart, zend) in reference coordinates | ||
| 7 | !> - Cartesian positions (xstart_cart, xend_cart) for visualization/heat loads | ||
| 8 | !> - Simulation config as global attributes | ||
| 9 | !> | ||
| 10 | !> Usage: | ||
| 11 | !> call write_results_netcdf(filename) | ||
| 12 | !> | ||
| 13 | !> The Cartesian positions decouple post-processing from coordinate | ||
| 14 | !> transformation machinery - Python can directly read (x,y,z) positions. | ||
| 15 | |||
| 16 | use, intrinsic :: iso_fortran_env, only: dp => real64, int8 | ||
| 17 | use netcdf | ||
| 18 | implicit none | ||
| 19 | |||
| 20 | private | ||
| 21 | public :: write_results_netcdf | ||
| 22 | |||
| 23 | contains | ||
| 24 | |||
| 25 | 145 | subroutine check_nc(status, location) | |
| 26 | !> Check NetCDF status and stop on error. | ||
| 27 | integer, intent(in) :: status | ||
| 28 | character(len=*), intent(in) :: location | ||
| 29 | |||
| 30 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 144 times.
|
144 | if (status /= nf90_noerr) then |
| 31 | ✗ | print *, 'NetCDF error at ', trim(location), ': ', & | |
| 32 | ✗ | trim(nf90_strerror(status)) | |
| 33 | ✗ | error stop 'NetCDF operation failed' | |
| 34 | end if | ||
| 35 | 144 | end subroutine check_nc | |
| 36 | |||
| 37 | 1 | subroutine write_results_netcdf(filename, wall_facet_count) | |
| 38 | !> Write all particle results to NetCDF file. | ||
| 39 | !> | ||
| 40 | !> Computes Cartesian positions from reference coordinates and writes | ||
| 41 | !> complete results including config attributes. | ||
| 42 | |||
| 43 | use params, only: ntestpart, trace_time, zstart, zend, times_lost, & | ||
| 44 | orbit_exit_code, boundary_event_radial_residual, & | ||
| 45 | boundary_event_time_width, trap_par, perp_inv, iclass, & | ||
| 46 | class_lost, & | ||
| 47 | facE_al, particle_energy_eV, & | ||
| 48 | particle_energy_eV_effective, n_d, n_e, v0, & | ||
| 49 | collision_model, nu_star_standard, & | ||
| 50 | lorentz_major_radius_cm, lorentz_iota, lorentz_nu, & | ||
| 51 | radial_electric_potential_slope, & | ||
| 52 | integmode, relerr, npoiper, npoiper2, & | ||
| 53 | boundary_event_fraction_tolerance, & | ||
| 54 | boundary_event_radial_tolerance, & | ||
| 55 | canonical_grid_nr, canonical_grid_ntheta, & | ||
| 56 | canonical_grid_nphi, canonical_ode_relerr, & | ||
| 57 | ntimstep, startmode, num_surf, sbeg, & | ||
| 58 | contr_pp, notrace_passing, & | ||
| 59 | isw_field_type, swcoll, deterministic, ran_seed, & | ||
| 60 | max_consecutive_warning_holds, & | ||
| 61 | netcdffile, field_input, coord_input, & | ||
| 62 | wall_input, wall_units, wall_hit, wall_hit_triangle, & | ||
| 63 | wall_hit_cart, & | ||
| 64 | wall_query_rho_lcfs, & | ||
| 65 | chart_boundary_kind, chart_boundary_kind_effective, & | ||
| 66 | wall_hit_normal_cart, wall_hit_cos_incidence, & | ||
| 67 | wall_hit_angle_rad | ||
| 68 | use new_vmec_stuff_mod, only: multharm | ||
| 69 | use chartmap_metadata, only: chartmap_metadata_t, read_chartmap_metadata | ||
| 70 | use libneo_coordinates, only: chartmap_coordinate_system_t | ||
| 71 | use reference_coordinates, only: ref_coords | ||
| 72 | use version, only: simple_version | ||
| 73 | use diag_counters, only: N_EVENT, diag_counters_total, event_name | ||
| 74 | use orbit_symplectic_base, only: symplectic_newton_warning_factor, & | ||
| 75 | symplectic_spectre_newton_warning_factor | ||
| 76 | |||
| 77 | character(len=*), intent(in) :: filename | ||
| 78 | integer, intent(in), optional :: wall_facet_count | ||
| 79 | |||
| 80 | integer :: event_id, ncid, status | ||
| 81 | 1 | character(:), allocatable :: diagnostic_name | |
| 82 | integer :: dim_particle, dim_xyz, dim_phase, dim_iclass | ||
| 83 | integer :: var_times_lost, var_orbit_exit_code, var_trap_par, var_perp_inv | ||
| 84 | integer :: var_boundary_event_radial_residual, var_boundary_event_time_width | ||
| 85 | integer :: var_zstart, var_zend, var_xstart_cart, var_xend_cart | ||
| 86 | integer :: var_iclass, var_class_lost, var_wall_hit, var_wall_hit_cart | ||
| 87 | integer :: var_wall_hit_normal_cart, var_wall_hit_cos_incidence | ||
| 88 | integer :: var_wall_hit_angle_rad, var_wall_hit_triangle | ||
| 89 | integer :: i | ||
| 90 | 1 | real(dp), allocatable :: xstart_cart(:, :), xend_cart(:, :) | |
| 91 | 1 | integer(int8), allocatable :: class_lost_i8(:) | |
| 92 | character(len=32) :: coord_type | ||
| 93 | character(len=16) :: cart_units | ||
| 94 | type(chartmap_metadata_t) :: meta | ||
| 95 | character(len=8) :: date_str | ||
| 96 | character(len=10) :: time_str | ||
| 97 | character(len=64) :: timestamp | ||
| 98 | |||
| 99 | ! Allocate working arrays | ||
| 100 |
5/10✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 time.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 time.
|
3 | allocate (xstart_cart(3, ntestpart)) |
| 101 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
|
2 | allocate (xend_cart(3, ntestpart)) |
| 102 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 time.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 time.
|
3 | allocate (class_lost_i8(ntestpart)) |
| 103 | |||
| 104 | ! Get current timestamp | ||
| 105 | 1 | call date_and_time(date=date_str, time=time_str) | |
| 106 | write (timestamp, '(A4,"-",A2,"-",A2,"T",A2,":",A2,":",A2)') & | ||
| 107 | 1 | date_str(1:4), date_str(5:6), date_str(7:8), & | |
| 108 | 2 | time_str(1:2), time_str(3:4), time_str(5:6) | |
| 109 | |||
| 110 | ! Compute Cartesian positions | ||
| 111 | 1 | call compute_cartesian_positions(xstart_cart, xend_cart) | |
| 112 | |||
| 113 | ! Convert logical to int8 | ||
| 114 |
2/2✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 time.
|
33 | do i = 1, ntestpart |
| 115 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
|
33 | if (class_lost(i)) then |
| 116 | ✗ | class_lost_i8(i) = 1_int8 | |
| 117 | else | ||
| 118 | 32 | class_lost_i8(i) = 0_int8 | |
| 119 | end if | ||
| 120 | end do | ||
| 121 | |||
| 122 | 1 | coord_type = 'vmec' | |
| 123 | 1 | cart_units = 'cm' | |
| 124 | select type (ref_coords) | ||
| 125 | class is (chartmap_coordinate_system_t) | ||
| 126 | ✗ | coord_type = 'chartmap' | |
| 127 | ✗ | call read_chartmap_metadata(coord_input, meta) | |
| 128 | ✗ | cart_units = meta%cart_units | |
| 129 | class default | ||
| 130 | continue | ||
| 131 | end select | ||
| 132 | |||
| 133 | ! Create file | ||
| 134 | 1 | status = nf90_create(filename, nf90_netcdf4, ncid) | |
| 135 | 1 | call check_nc(status, 'nf90_create') | |
| 136 | |||
| 137 | ! Define dimensions | ||
| 138 | 1 | status = nf90_def_dim(ncid, 'particle', ntestpart, dim_particle) | |
| 139 | 1 | call check_nc(status, 'def_dim particle') | |
| 140 | |||
| 141 | 1 | status = nf90_def_dim(ncid, 'xyz', 3, dim_xyz) | |
| 142 | 1 | call check_nc(status, 'def_dim xyz') | |
| 143 | |||
| 144 | 1 | status = nf90_def_dim(ncid, 'phase', 5, dim_phase) | |
| 145 | 1 | call check_nc(status, 'def_dim phase') | |
| 146 | |||
| 147 | 1 | status = nf90_def_dim(ncid, 'iclass_dim', 3, dim_iclass) | |
| 148 | 1 | call check_nc(status, 'def_dim iclass_dim') | |
| 149 | |||
| 150 | ! Define variables with full attribute error checking | ||
| 151 | status = nf90_def_var(ncid, 'times_lost', nf90_double, & | ||
| 152 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_times_lost) |
| 153 | 1 | call check_nc(status, 'def_var times_lost') | |
| 154 | 1 | status = nf90_put_att(ncid, var_times_lost, 'units', 's') | |
| 155 | 1 | call check_nc(status, 'put_att times_lost units') | |
| 156 | status = nf90_put_att(ncid, var_times_lost, 'description', & | ||
| 157 | 'Physical loss time (-1 if skipped; trace_time if confined; ' // & | ||
| 158 | 1 | 'NaN after fatal numerical failure)') | |
| 159 | 1 | call check_nc(status, 'put_att times_lost description') | |
| 160 | |||
| 161 | status = nf90_def_var(ncid, 'orbit_exit_code', nf90_int, & | ||
| 162 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_orbit_exit_code) |
| 163 | 1 | call check_nc(status, 'def_var orbit_exit_code') | |
| 164 | status = nf90_put_att(ncid, var_orbit_exit_code, 'description', & | ||
| 165 | '0=completed, 1=LCFS, 2=wall, 3=skipped, ' // & | ||
| 166 | 1 | '4=numerically confined in core, 101-105=numerical failure') | |
| 167 | 1 | call check_nc(status, 'put_att orbit_exit_code description') | |
| 168 | |||
| 169 | status = nf90_def_var(ncid, 'boundary_event_radial_residual', nf90_double, & | ||
| 170 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_boundary_event_radial_residual) |
| 171 | 1 | call check_nc(status, 'def_var boundary_event_radial_residual') | |
| 172 | status = nf90_put_att(ncid, var_boundary_event_radial_residual, 'units', & | ||
| 173 | 1 | 'normalized radial coordinate') | |
| 174 | 1 | call check_nc(status, 'put_att boundary_event_radial_residual units') | |
| 175 | |||
| 176 | status = nf90_def_var(ncid, 'boundary_event_time_width', nf90_double, & | ||
| 177 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_boundary_event_time_width) |
| 178 | 1 | call check_nc(status, 'def_var boundary_event_time_width') | |
| 179 | 1 | status = nf90_put_att(ncid, var_boundary_event_time_width, 'units', 's') | |
| 180 | 1 | call check_nc(status, 'put_att boundary_event_time_width units') | |
| 181 | |||
| 182 | status = nf90_def_var(ncid, 'trap_par', nf90_double, & | ||
| 183 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_trap_par) |
| 184 | 1 | call check_nc(status, 'def_var trap_par') | |
| 185 | 1 | status = nf90_put_att(ncid, var_trap_par, 'units', 'dimensionless') | |
| 186 | 1 | call check_nc(status, 'put_att trap_par units') | |
| 187 | status = nf90_put_att(ncid, var_trap_par, 'description', & | ||
| 188 | 'Trapping parameter (mu*B_max/E - 1)*B_min/' // & | ||
| 189 | '(B_max - B_min) on the start surface: ' // & | ||
| 190 | '<0 passing, 0 trapped-passing boundary, ' // & | ||
| 191 | 1 | '1 reflecting at B_min') | |
| 192 | 1 | call check_nc(status, 'put_att trap_par description') | |
| 193 | |||
| 194 | status = nf90_def_var(ncid, 'perp_inv', nf90_double, & | ||
| 195 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_perp_inv) |
| 196 | 1 | call check_nc(status, 'def_var perp_inv') | |
| 197 | 1 | status = nf90_put_att(ncid, var_perp_inv, 'units', 'cm^2*G') | |
| 198 | 1 | call check_nc(status, 'put_att perp_inv units') | |
| 199 | status = nf90_put_att(ncid, var_perp_inv, 'description', & | ||
| 200 | 1 | 'Perpendicular adiabatic invariant J_perp = mu*B') | |
| 201 | 1 | call check_nc(status, 'put_att perp_inv description') | |
| 202 | |||
| 203 | status = nf90_def_var(ncid, 'zstart', nf90_double, & | ||
| 204 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | [dim_phase, dim_particle], var_zstart) |
| 205 | 1 | call check_nc(status, 'def_var zstart') | |
| 206 | status = nf90_put_att(ncid, var_zstart, 'description', & | ||
| 207 | 1 | 'Initial phase space in reference coordinates') | |
| 208 | 1 | call check_nc(status, 'put_att zstart description') | |
| 209 | status = nf90_put_att(ncid, var_zstart, 'components', & | ||
| 210 | 1 | 's (flux), theta (rad), zeta (rad), p (v/v0), v_par/v') | |
| 211 | 1 | call check_nc(status, 'put_att zstart components') | |
| 212 | status = nf90_put_att(ncid, var_zstart, 'coordinate_note', & | ||
| 213 | 1 | 'See coordinate_type attribute for VMEC vs chartmap') | |
| 214 | 1 | call check_nc(status, 'put_att zstart coordinate_note') | |
| 215 | |||
| 216 | status = nf90_def_var(ncid, 'zend', nf90_double, & | ||
| 217 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | [dim_phase, dim_particle], var_zend) |
| 218 | 1 | call check_nc(status, 'def_var zend') | |
| 219 | status = nf90_put_att(ncid, var_zend, 'description', & | ||
| 220 | 1 | 'Final phase space in ref coords (0 if not traced)') | |
| 221 | 1 | call check_nc(status, 'put_att zend description') | |
| 222 | status = nf90_put_att(ncid, var_zend, 'components', & | ||
| 223 | 1 | 's (flux), theta (rad), zeta (rad), p (v/v0), v_par/v') | |
| 224 | 1 | call check_nc(status, 'put_att zend components') | |
| 225 | |||
| 226 | status = nf90_def_var(ncid, 'xstart_cart', nf90_double, & | ||
| 227 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | [dim_xyz, dim_particle], var_xstart_cart) |
| 228 | 1 | call check_nc(status, 'def_var xstart_cart') | |
| 229 | 1 | status = nf90_put_att(ncid, var_xstart_cart, 'units', trim(cart_units)) | |
| 230 | 1 | call check_nc(status, 'put_att xstart_cart units') | |
| 231 | status = nf90_put_att(ncid, var_xstart_cart, 'description', & | ||
| 232 | 1 | 'Initial Cartesian position (x, y, z)') | |
| 233 | 1 | call check_nc(status, 'put_att xstart_cart description') | |
| 234 | |||
| 235 | status = nf90_def_var(ncid, 'xend_cart', nf90_double, & | ||
| 236 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | [dim_xyz, dim_particle], var_xend_cart) |
| 237 | 1 | call check_nc(status, 'def_var xend_cart') | |
| 238 | 1 | status = nf90_put_att(ncid, var_xend_cart, 'units', trim(cart_units)) | |
| 239 | 1 | call check_nc(status, 'put_att xend_cart units') | |
| 240 | status = nf90_put_att(ncid, var_xend_cart, 'description', & | ||
| 241 | 1 | 'Final Cartesian position (= xstart if not traced)') | |
| 242 | 1 | call check_nc(status, 'put_att xend_cart description') | |
| 243 | |||
| 244 | status = nf90_def_var(ncid, 'iclass', nf90_int, & | ||
| 245 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | [dim_iclass, dim_particle], var_iclass) |
| 246 | 1 | call check_nc(status, 'def_var iclass') | |
| 247 | status = nf90_put_att(ncid, var_iclass, 'description', & | ||
| 248 | 'Classification: [J_par, topology, fractal] '// & | ||
| 249 | 1 | '(0=unclassified)') | |
| 250 | 1 | call check_nc(status, 'put_att iclass description') | |
| 251 | |||
| 252 | status = nf90_def_var(ncid, 'class_lost', nf90_byte, & | ||
| 253 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_class_lost) |
| 254 | 1 | call check_nc(status, 'def_var class_lost') | |
| 255 | status = nf90_put_att(ncid, var_class_lost, 'description', & | ||
| 256 | 1 | 'Lost flag (1=lost, 0=confined)') | |
| 257 | 1 | call check_nc(status, 'put_att class_lost description') | |
| 258 | |||
| 259 | status = nf90_def_var(ncid, 'wall_hit', nf90_byte, & | ||
| 260 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_wall_hit) |
| 261 | 1 | call check_nc(status, 'def_var wall_hit') | |
| 262 | status = nf90_put_att(ncid, var_wall_hit, 'description', & | ||
| 263 | 1 | 'Wall hit flag (1=hit wall_input STL, 0=no hit)') | |
| 264 | 1 | call check_nc(status, 'put_att wall_hit description') | |
| 265 | |||
| 266 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (len_trim(wall_input) > 0) then |
| 267 | status = nf90_def_var(ncid, 'wall_hit_triangle', nf90_int, & | ||
| 268 | ✗ | [dim_particle], var_wall_hit_triangle) | |
| 269 | ✗ | call check_nc(status, 'def_var wall_hit_triangle') | |
| 270 | status = nf90_put_att(ncid, var_wall_hit_triangle, 'description', & | ||
| 271 | ✗ | 'STL facet index in file order; -1 means no wall hit') | |
| 272 | ✗ | call check_nc(status, 'put_att wall_hit_triangle description') | |
| 273 | ✗ | status = nf90_put_att(ncid, var_wall_hit_triangle, 'start_index', 0) | |
| 274 | ✗ | call check_nc(status, 'put_att wall_hit_triangle start_index') | |
| 275 | status = nf90_put_att(ncid, var_wall_hit_triangle, 'stl_file', & | ||
| 276 | ✗ | trim(wall_input)) | |
| 277 | ✗ | call check_nc(status, 'put_att wall_hit_triangle stl_file') | |
| 278 | ✗ | if (.not. present(wall_facet_count)) & | |
| 279 | ✗ | error stop 'STL output requires the initialized wall facet count' | |
| 280 | status = nf90_put_att(ncid, var_wall_hit_triangle, & | ||
| 281 | ✗ | 'stl_facet_count', wall_facet_count) | |
| 282 | ✗ | call check_nc(status, 'put_att wall_hit_triangle stl_facet_count') | |
| 283 | end if | ||
| 284 | |||
| 285 | status = nf90_def_var(ncid, 'wall_hit_cart', nf90_double, & | ||
| 286 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | [dim_xyz, dim_particle], var_wall_hit_cart) |
| 287 | 1 | call check_nc(status, 'def_var wall_hit_cart') | |
| 288 | 1 | status = nf90_put_att(ncid, var_wall_hit_cart, 'units', trim(cart_units)) | |
| 289 | 1 | call check_nc(status, 'put_att wall_hit_cart units') | |
| 290 | status = nf90_put_att(ncid, var_wall_hit_cart, 'description', & | ||
| 291 | 1 | 'First-hit Cartesian wall intersection point (x, y, z)') | |
| 292 | 1 | call check_nc(status, 'put_att wall_hit_cart description') | |
| 293 | |||
| 294 | status = nf90_def_var(ncid, 'wall_hit_normal_cart', nf90_double, & | ||
| 295 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | [dim_xyz, dim_particle], var_wall_hit_normal_cart) |
| 296 | 1 | call check_nc(status, 'def_var wall_hit_normal_cart') | |
| 297 | 1 | status = nf90_put_att(ncid, var_wall_hit_normal_cart, 'units', '1') | |
| 298 | 1 | call check_nc(status, 'put_att wall_hit_normal_cart units') | |
| 299 | status = nf90_put_att(ncid, var_wall_hit_normal_cart, 'description', & | ||
| 300 | 1 | 'Unit surface normal at wall hit (x, y, z)') | |
| 301 | 1 | call check_nc(status, 'put_att wall_hit_normal_cart description') | |
| 302 | |||
| 303 | status = nf90_def_var(ncid, 'wall_hit_cos_incidence', nf90_double, & | ||
| 304 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_wall_hit_cos_incidence) |
| 305 | 1 | call check_nc(status, 'def_var wall_hit_cos_incidence') | |
| 306 | 1 | status = nf90_put_att(ncid, var_wall_hit_cos_incidence, 'units', '1') | |
| 307 | 1 | call check_nc(status, 'put_att wall_hit_cos_incidence units') | |
| 308 | status = nf90_put_att(ncid, var_wall_hit_cos_incidence, 'description', & | ||
| 309 | 1 | 'Absolute cosine of incidence angle at wall hit') | |
| 310 | 1 | call check_nc(status, 'put_att wall_hit_cos_incidence description') | |
| 311 | |||
| 312 | status = nf90_def_var(ncid, 'wall_hit_angle_rad', nf90_double, & | ||
| 313 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | [dim_particle], var_wall_hit_angle_rad) |
| 314 | 1 | call check_nc(status, 'def_var wall_hit_angle_rad') | |
| 315 | 1 | status = nf90_put_att(ncid, var_wall_hit_angle_rad, 'units', 'rad') | |
| 316 | 1 | call check_nc(status, 'put_att wall_hit_angle_rad units') | |
| 317 | status = nf90_put_att(ncid, var_wall_hit_angle_rad, 'description', & | ||
| 318 | 1 | 'Incidence angle at wall hit (0 = normal incidence)') | |
| 319 | 1 | call check_nc(status, 'put_att wall_hit_angle_rad description') | |
| 320 | |||
| 321 | ! Global attributes - simulation config | ||
| 322 | status = nf90_put_att(ncid, nf90_global, 'title', & | ||
| 323 | 1 | 'SIMPLE particle tracing results') | |
| 324 | 1 | call check_nc(status, 'put_att title') | |
| 325 | status = nf90_put_att(ncid, nf90_global, 'simple_version', & | ||
| 326 | 1 | trim(simple_version)) | |
| 327 | 1 | call check_nc(status, 'put_att simple_version') | |
| 328 | 1 | status = nf90_put_att(ncid, nf90_global, 'cart_units', trim(cart_units)) | |
| 329 | 1 | call check_nc(status, 'put_att cart_units') | |
| 330 | 1 | status = nf90_put_att(ncid, nf90_global, 'created', trim(timestamp)) | |
| 331 | 1 | call check_nc(status, 'put_att created') | |
| 332 | 1 | status = nf90_put_att(ncid, nf90_global, 'ntestpart', ntestpart) | |
| 333 | 1 | call check_nc(status, 'put_att ntestpart') | |
| 334 | 1 | status = nf90_put_att(ncid, nf90_global, 'trace_time', trace_time) | |
| 335 | 1 | call check_nc(status, 'put_att trace_time') | |
| 336 | 1 | status = nf90_put_att(ncid, nf90_global, 'facE_al', facE_al) | |
| 337 | 1 | call check_nc(status, 'put_att facE_al') | |
| 338 | status = nf90_put_att(ncid, nf90_global, 'particle_energy_eV_input', & | ||
| 339 | 1 | particle_energy_eV) | |
| 340 | 1 | call check_nc(status, 'put_att particle_energy_eV_input') | |
| 341 | status = nf90_put_att(ncid, nf90_global, 'particle_energy_eV', & | ||
| 342 | 1 | particle_energy_eV_effective) | |
| 343 | 1 | call check_nc(status, 'put_att particle_energy_eV') | |
| 344 | 1 | status = nf90_put_att(ncid, nf90_global, 'particle_mass_number', n_d) | |
| 345 | 1 | call check_nc(status, 'put_att particle_mass_number') | |
| 346 | 1 | status = nf90_put_att(ncid, nf90_global, 'particle_charge_number', n_e) | |
| 347 | 1 | call check_nc(status, 'put_att particle_charge_number') | |
| 348 | 1 | status = nf90_put_att(ncid, nf90_global, 'v0', v0) | |
| 349 | 1 | call check_nc(status, 'put_att v0') | |
| 350 | status = nf90_put_att(ncid, nf90_global, 'collision_model', & | ||
| 351 | 1 | trim(collision_model)) | |
| 352 | 1 | call check_nc(status, 'put_att collision_model') | |
| 353 | status = nf90_put_att(ncid, nf90_global, 'nu_star_standard', & | ||
| 354 | 1 | nu_star_standard) | |
| 355 | 1 | call check_nc(status, 'put_att nu_star_standard') | |
| 356 | status = nf90_put_att(ncid, nf90_global, 'lorentz_major_radius_cm', & | ||
| 357 | 1 | lorentz_major_radius_cm) | |
| 358 | 1 | call check_nc(status, 'put_att lorentz_major_radius_cm') | |
| 359 | 1 | status = nf90_put_att(ncid, nf90_global, 'lorentz_iota', lorentz_iota) | |
| 360 | 1 | call check_nc(status, 'put_att lorentz_iota') | |
| 361 | 1 | status = nf90_put_att(ncid, nf90_global, 'lorentz_nu_s-1', lorentz_nu) | |
| 362 | 1 | call check_nc(status, 'put_att lorentz_nu_s-1') | |
| 363 | status = nf90_put_att(ncid, nf90_global, & | ||
| 364 | 1 | 'radial_electric_potential_slope', radial_electric_potential_slope) | |
| 365 | 1 | call check_nc(status, 'put_att radial_electric_potential_slope') | |
| 366 | 1 | status = nf90_put_att(ncid, nf90_global, 'integmode', integmode) | |
| 367 | 1 | call check_nc(status, 'put_att integmode') | |
| 368 | 1 | status = nf90_put_att(ncid, nf90_global, 'relerr', relerr) | |
| 369 | 1 | call check_nc(status, 'put_att relerr') | |
| 370 | status = nf90_put_att(ncid, nf90_global, & | ||
| 371 | 'symplectic_newton_warning_factor', & | ||
| 372 | 1 | symplectic_newton_warning_factor) | |
| 373 | 1 | call check_nc(status, 'put_att symplectic_newton_warning_factor') | |
| 374 | status = nf90_put_att(ncid, nf90_global, & | ||
| 375 | 'symplectic_spectre_newton_warning_factor', & | ||
| 376 | 1 | symplectic_spectre_newton_warning_factor) | |
| 377 | 1 | call check_nc(status, 'put_att symplectic_spectre_newton_warning_factor') | |
| 378 | status = nf90_put_att(ncid, nf90_global, & | ||
| 379 | 'boundary_event_fraction_tolerance', & | ||
| 380 | 1 | boundary_event_fraction_tolerance) | |
| 381 | 1 | call check_nc(status, 'put_att boundary_event_fraction_tolerance') | |
| 382 | status = nf90_put_att(ncid, nf90_global, & | ||
| 383 | 'boundary_event_radial_tolerance', & | ||
| 384 | 1 | boundary_event_radial_tolerance) | |
| 385 | 1 | call check_nc(status, 'put_att boundary_event_radial_tolerance') | |
| 386 | status = nf90_put_att(ncid, nf90_global, 'canonical_grid_nr', & | ||
| 387 | 1 | canonical_grid_nr) | |
| 388 | 1 | call check_nc(status, 'put_att canonical_grid_nr') | |
| 389 | status = nf90_put_att(ncid, nf90_global, 'canonical_grid_ntheta', & | ||
| 390 | 1 | canonical_grid_ntheta) | |
| 391 | 1 | call check_nc(status, 'put_att canonical_grid_ntheta') | |
| 392 | status = nf90_put_att(ncid, nf90_global, 'canonical_grid_nphi', & | ||
| 393 | 1 | canonical_grid_nphi) | |
| 394 | 1 | call check_nc(status, 'put_att canonical_grid_nphi') | |
| 395 | status = nf90_put_att(ncid, nf90_global, 'canonical_ode_relerr', & | ||
| 396 | 1 | canonical_ode_relerr) | |
| 397 | 1 | call check_nc(status, 'put_att canonical_ode_relerr') | |
| 398 | 1 | status = nf90_put_att(ncid, nf90_global, 'multharm', multharm) | |
| 399 | 1 | call check_nc(status, 'put_att multharm') | |
| 400 | 1 | status = nf90_put_att(ncid, nf90_global, 'npoiper', npoiper) | |
| 401 | 1 | call check_nc(status, 'put_att npoiper') | |
| 402 | 1 | status = nf90_put_att(ncid, nf90_global, 'npoiper2', npoiper2) | |
| 403 | 1 | call check_nc(status, 'put_att npoiper2') | |
| 404 | 1 | status = nf90_put_att(ncid, nf90_global, 'ntimstep', ntimstep) | |
| 405 | 1 | call check_nc(status, 'put_att ntimstep') | |
| 406 | 1 | status = nf90_put_att(ncid, nf90_global, 'startmode', startmode) | |
| 407 | 1 | call check_nc(status, 'put_att startmode') | |
| 408 | 1 | status = nf90_put_att(ncid, nf90_global, 'num_surf', num_surf) | |
| 409 | 1 | call check_nc(status, 'put_att num_surf') | |
| 410 | 1 | status = nf90_put_att(ncid, nf90_global, 'sbeg', sbeg) | |
| 411 | 1 | call check_nc(status, 'put_att sbeg') | |
| 412 | 1 | status = nf90_put_att(ncid, nf90_global, 'contr_pp', contr_pp) | |
| 413 | 1 | call check_nc(status, 'put_att contr_pp') | |
| 414 | status = nf90_put_att(ncid, nf90_global, 'notrace_passing', & | ||
| 415 | 1 | notrace_passing) | |
| 416 | 1 | call check_nc(status, 'put_att notrace_passing') | |
| 417 | 1 | status = nf90_put_att(ncid, nf90_global, 'isw_field_type', isw_field_type) | |
| 418 | 1 | call check_nc(status, 'put_att isw_field_type') | |
| 419 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
2 | status = nf90_put_att(ncid, nf90_global, 'swcoll', merge(1, 0, swcoll)) |
| 420 | 1 | call check_nc(status, 'put_att swcoll') | |
| 421 | status = nf90_put_att(ncid, nf90_global, 'deterministic', & | ||
| 422 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | merge(1, 0, deterministic)) |
| 423 | 1 | call check_nc(status, 'put_att deterministic') | |
| 424 | 1 | status = nf90_put_att(ncid, nf90_global, 'ran_seed', ran_seed) | |
| 425 | 1 | call check_nc(status, 'put_att ran_seed') | |
| 426 | status = nf90_put_att(ncid, nf90_global, & | ||
| 427 | 1 | 'max_consecutive_warning_holds', max_consecutive_warning_holds) | |
| 428 | 1 | call check_nc(status, 'put_att max_consecutive_warning_holds') | |
| 429 | 1 | status = nf90_put_att(ncid, nf90_global, 'netcdffile', trim(netcdffile)) | |
| 430 | 1 | call check_nc(status, 'put_att netcdffile') | |
| 431 | 1 | status = nf90_put_att(ncid, nf90_global, 'field_input', trim(field_input)) | |
| 432 | 1 | call check_nc(status, 'put_att field_input') | |
| 433 | 1 | status = nf90_put_att(ncid, nf90_global, 'coord_input', trim(coord_input)) | |
| 434 | 1 | call check_nc(status, 'put_att coord_input') | |
| 435 | 1 | status = nf90_put_att(ncid, nf90_global, 'wall_input', trim(wall_input)) | |
| 436 | 1 | call check_nc(status, 'put_att wall_input') | |
| 437 | 1 | status = nf90_put_att(ncid, nf90_global, 'wall_units', trim(wall_units)) | |
| 438 | 1 | call check_nc(status, 'put_att wall_units') | |
| 439 | status = nf90_put_att(ncid, nf90_global, 'wall_query_rho_lcfs', & | ||
| 440 | 1 | wall_query_rho_lcfs) | |
| 441 | 1 | call check_nc(status, 'put_att wall_query_rho_lcfs') | |
| 442 | status = nf90_put_att(ncid, nf90_global, 'chart_boundary_kind', & | ||
| 443 | 1 | trim(chart_boundary_kind)) | |
| 444 | 1 | call check_nc(status, 'put_att chart_boundary_kind') | |
| 445 | status = nf90_put_att(ncid, nf90_global, & | ||
| 446 | 'chart_boundary_kind_effective', & | ||
| 447 | 1 | trim(chart_boundary_kind_effective)) | |
| 448 | 1 | call check_nc(status, 'put_att chart_boundary_kind_effective') | |
| 449 | 1 | status = nf90_put_att(ncid, nf90_global, 'coordinate_type', trim(coord_type)) | |
| 450 | 1 | call check_nc(status, 'put_att coordinate_type') | |
| 451 |
2/2✓ Branch 0 taken 23 times.
✓ Branch 1 taken 1 time.
|
24 | do event_id = 1, N_EVENT |
| 452 |
5/6✓ Branch 2 taken 1 time.
✓ Branch 3 taken 22 times.
✓ Branch 4 taken 20 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 23 times.
✗ Branch 7 not taken.
|
23 | diagnostic_name = 'diagnostic_'//event_name(event_id) |
| 453 | status = nf90_put_att(ncid, nf90_global, diagnostic_name, & | ||
| 454 | 23 | diag_counters_total(event_id)) | |
| 455 | 24 | call check_nc(status, 'put_att '//diagnostic_name) | |
| 456 | end do | ||
| 457 | |||
| 458 | ! End define mode | ||
| 459 | 1 | status = nf90_enddef(ncid) | |
| 460 | 1 | call check_nc(status, 'enddef') | |
| 461 | |||
| 462 | ! Write data | ||
| 463 | 1 | status = nf90_put_var(ncid, var_times_lost, times_lost) | |
| 464 | 1 | call check_nc(status, 'put_var times_lost') | |
| 465 | |||
| 466 | 1 | status = nf90_put_var(ncid, var_orbit_exit_code, orbit_exit_code) | |
| 467 | 1 | call check_nc(status, 'put_var orbit_exit_code') | |
| 468 | |||
| 469 | status = nf90_put_var(ncid, var_boundary_event_radial_residual, & | ||
| 470 | 1 | boundary_event_radial_residual) | |
| 471 | 1 | call check_nc(status, 'put_var boundary_event_radial_residual') | |
| 472 | |||
| 473 | status = nf90_put_var(ncid, var_boundary_event_time_width, & | ||
| 474 | 1 | boundary_event_time_width) | |
| 475 | 1 | call check_nc(status, 'put_var boundary_event_time_width') | |
| 476 | |||
| 477 | 1 | status = nf90_put_var(ncid, var_trap_par, trap_par) | |
| 478 | 1 | call check_nc(status, 'put_var trap_par') | |
| 479 | |||
| 480 | 1 | status = nf90_put_var(ncid, var_perp_inv, perp_inv) | |
| 481 | 1 | call check_nc(status, 'put_var perp_inv') | |
| 482 | |||
| 483 | 1 | status = nf90_put_var(ncid, var_zstart, zstart) | |
| 484 | 1 | call check_nc(status, 'put_var zstart') | |
| 485 | |||
| 486 | 1 | status = nf90_put_var(ncid, var_zend, zend) | |
| 487 | 1 | call check_nc(status, 'put_var zend') | |
| 488 | |||
| 489 | 1 | status = nf90_put_var(ncid, var_xstart_cart, xstart_cart) | |
| 490 | 1 | call check_nc(status, 'put_var xstart_cart') | |
| 491 | |||
| 492 | 1 | status = nf90_put_var(ncid, var_xend_cart, xend_cart) | |
| 493 | 1 | call check_nc(status, 'put_var xend_cart') | |
| 494 | |||
| 495 | 1 | status = nf90_put_var(ncid, var_iclass, iclass) | |
| 496 | 1 | call check_nc(status, 'put_var iclass') | |
| 497 | |||
| 498 | 1 | status = nf90_put_var(ncid, var_class_lost, class_lost_i8) | |
| 499 | 1 | call check_nc(status, 'put_var class_lost') | |
| 500 | |||
| 501 | 1 | status = nf90_put_var(ncid, var_wall_hit, wall_hit) | |
| 502 | 1 | call check_nc(status, 'put_var wall_hit') | |
| 503 | |||
| 504 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (len_trim(wall_input) > 0) then |
| 505 | ✗ | status = nf90_put_var(ncid, var_wall_hit_triangle, wall_hit_triangle) | |
| 506 | ✗ | call check_nc(status, 'put_var wall_hit_triangle') | |
| 507 | end if | ||
| 508 | |||
| 509 | 1 | status = nf90_put_var(ncid, var_wall_hit_cart, wall_hit_cart) | |
| 510 | 1 | call check_nc(status, 'put_var wall_hit_cart') | |
| 511 | |||
| 512 | 1 | status = nf90_put_var(ncid, var_wall_hit_normal_cart, wall_hit_normal_cart) | |
| 513 | 1 | call check_nc(status, 'put_var wall_hit_normal_cart') | |
| 514 | |||
| 515 | 1 | status = nf90_put_var(ncid, var_wall_hit_cos_incidence, wall_hit_cos_incidence) | |
| 516 | 1 | call check_nc(status, 'put_var wall_hit_cos_incidence') | |
| 517 | |||
| 518 | 1 | status = nf90_put_var(ncid, var_wall_hit_angle_rad, wall_hit_angle_rad) | |
| 519 | 1 | call check_nc(status, 'put_var wall_hit_angle_rad') | |
| 520 | |||
| 521 | ! Close file | ||
| 522 | 1 | status = nf90_close(ncid) | |
| 523 | 1 | call check_nc(status, 'close') | |
| 524 | |||
| 525 | 1 | print *, 'INFO: Results written to ', trim(filename) | |
| 526 | |||
| 527 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 time.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✗ Branch 9 not taken.
|
2 | end subroutine write_results_netcdf |
| 528 | |||
| 529 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
|
1 | subroutine compute_cartesian_positions(xstart_cart, xend_cart) |
| 530 | !> Convert reference coordinates to Cartesian for all particles. | ||
| 531 | !> | ||
| 532 | !> For particles that were not traced (zend = 0), xend_cart is set to zero. | ||
| 533 | !> | ||
| 534 | !> For STL wall hits (wall_hit == 1), xend_cart is set directly from | ||
| 535 | !> wall_hit_cart since that is the authoritative CGAL intersection point. | ||
| 536 | !> This avoids round-trip coordinate conversion errors when from_cart | ||
| 537 | !> fails in ill-conditioned chartmap regions. | ||
| 538 | |||
| 539 | 1 | use params, only: ntestpart, zstart, zend, wall_hit, wall_hit_cart | |
| 540 | use reference_coordinates, only: ref_coords | ||
| 541 | |||
| 542 | real(dp), intent(out) :: xstart_cart(:, :) | ||
| 543 | real(dp), intent(out) :: xend_cart(:, :) | ||
| 544 | integer :: i | ||
| 545 | logical :: zend_is_zero | ||
| 546 | |||
| 547 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if (.not. allocated(ref_coords)) then |
| 548 | ✗ | print *, 'ERROR: ref_coords not allocated - cannot compute Cartesian' | |
| 549 | ✗ | print *, 'This indicates init_reference_coordinates was not called' | |
| 550 | ✗ | error stop 'NetCDF results output requires initialized coordinate system' | |
| 551 | end if | ||
| 552 | |||
| 553 |
2/2✓ Branch 0 taken 32 times.
✓ Branch 1 taken 1 time.
|
33 | do i = 1, ntestpart |
| 554 |
2/8✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 32 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
32 | call ref_coords%evaluate_cart(zstart(1:3, i), xstart_cart(:, i)) |
| 555 | |||
| 556 | ! For STL wall hits, use the authoritative CGAL intersection point | ||
| 557 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
|
33 | if (wall_hit(i) == 1) then |
| 558 | ✗ | xend_cart(:, i) = wall_hit_cart(:, i) | |
| 559 | else | ||
| 560 | ! Check if particle was traced: zend is set to exactly 0 for untraced | ||
| 561 | ! particles (see simple_main.f90:493 and classification.f90:94) | ||
| 562 |
2/4✓ Branch 0 taken 32 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 32 times.
|
32 | zend_is_zero = all(zend(1:3, i) == 0.0_dp) |
| 563 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
|
32 | if (zend_is_zero) then |
| 564 | ! Particle not traced - keep xend_cart as zero | ||
| 565 | ✗ | xend_cart(:, i) = 0.0_dp | |
| 566 | else | ||
| 567 |
2/8✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 32 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
32 | call ref_coords%evaluate_cart(zend(1:3, i), xend_cart(:, i)) |
| 568 | end if | ||
| 569 | end if | ||
| 570 | end do | ||
| 571 | |||
| 572 | 1 | end subroutine compute_cartesian_positions | |
| 573 | |||
| 574 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
|
1 | end module netcdf_results_output |
| 575 |