diff --git a/examples/classification/Makefile b/examples/classification/Makefile index 7ad84e15..326293d3 100644 --- a/examples/classification/Makefile +++ b/examples/classification/Makefile @@ -1,24 +1,31 @@ -# Example for plotting classification over radius s and perpendicular invariant Jperp - -FC = gfortran -FFLAGS = -g -fbacktrace +# Classification example reproducing JPP 2020 paper results +# (Albert, Kasilov, Kernbichler, J. Plasma Phys. 86, 815860201) +# +# Configurations: QI (Drevlak 2014), QH (Drevlak 2018), QA (Henneberg 2019) +# +# Targets: +# make run - run all 10 cases under /tmp/simple_classification +# make plot - generate all paper-style plots from results +# make clean - remove local output files and /tmp results + +SIMPLE_X = ../../build/simple.x +BASE_DIR = /tmp/simple_classification all: plot -plot: prompt1.dat regular1.dat stochastic1.dat prompt2.dat regular2.dat stochastic2.dat plot_classification.py - python plot_classification.py - -prompt1.dat regular1.dat stochastic1.dat prompt2.dat regular2.dat stochastic2.dat: class_parts.dat postprocess_class.x - ./postprocess_class.x +run: $(SIMPLE_X) + bash run_all.sh -class_parts.dat: simple.in wout_LandremanPaul2021_QH_reactorScale_lowres_reference.nc - ../../build/simple.x +plot: + source ../../.venv/bin/activate && python plot_paper_results.py $(BASE_DIR) -wout_LandremanPaul2021_QH_reactorScale_lowres_reference.nc: - wget https://github.com/hiddenSymmetries/simsopt/raw/master/tests/test_files/wout_LandremanPaul2021_QH_reactorScale_lowres_reference.nc +# Legacy volume classification (runs in-place with simple.in) +plot-local: class_parts.dat bminmax.dat + source ../../.venv/bin/activate && python plot_classification.py -postprocess_class.x: postprocess_class.f90 - $(FC) $(FFLAGS) -o postprocess_class.x postprocess_class.f90 +class_parts.dat bminmax.dat: simple.in + $(SIMPLE_X) clean: - rm -f *.x *.dat fort.* *.pdf + rm -f *.dat fort.* *.pdf + rm -rf $(BASE_DIR) diff --git a/examples/classification/README.md b/examples/classification/README.md deleted file mode 100644 index 591dfddd..00000000 --- a/examples/classification/README.md +++ /dev/null @@ -1,26 +0,0 @@ -class_parts.dat columns: - -1: particle index -2: starting radial coordinate s -3: proportional to perpendicular adiabatic invariant (magnetic moment) -4: result of j_parallel orbitclassifier -5: result of topological (ideal) orbit classifier -6: fractal dimension classification (regular=1, chaotic=2) - -Prompt losses: 0 -Regular / ideal: 1 -Chaotic / non-ideal: 2 - -according to classification.f90, check_orbit_type.f90 and simple_main.f90 - -tcut is used only for fractal dimension classification. Earlier it was recommended to use at t=0.1s. The new classifiers (topological, and j_parallel) work at t=0.015s. - -Recipe: - -0) Minimize fraction of chaotic orbits. (10-100 times faster than usual run) -1) Look at s=0.25 and s=0.5. See at which J_perp particles are chaotic on s=0.25 and optimize such that at these J_perp values at s=0.5 there are regular regions. (10-100 times faster than usual run) -3) Optimize for loss fraction, tracing only chaotic orbits to the end. (only 1-10 times faster than usual run) - -TODO: - -Enable mode 3) also for new classifiers. \ No newline at end of file diff --git a/examples/classification/configs/qa_fig8.in b/examples/classification/configs/qa_fig8.in new file mode 100644 index 00000000..87a2f516 --- /dev/null +++ b/examples/classification/configs/qa_fig8.in @@ -0,0 +1,15 @@ +&config +notrace_passing = 0 +ntestpart = 10000 +num_surf = 1 +sbeg = 0.6 +npoiper2 = 128 +netcdffile = 'wout_henneberg_qa.nc' +isw_field_type = 2 +trace_time = 1.0d0 +tcut = 1d-1 +class_plot = .True. +cut_in_per = 0.5d0 +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qa_s03.in b/examples/classification/configs/qa_s03.in new file mode 100644 index 00000000..30738a04 --- /dev/null +++ b/examples/classification/configs/qa_s03.in @@ -0,0 +1,14 @@ +&config +notrace_passing = 1 +ntestpart = 100000 +num_surf = 1 +sbeg = 0.3 +npoiper2 = 128 +netcdffile = 'wout_henneberg_qa.nc' +isw_field_type = 2 +trace_time = 1.0d0 +tcut = 1d-1 +class_plot = .False. +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qa_s06.in b/examples/classification/configs/qa_s06.in new file mode 100644 index 00000000..737e52e7 --- /dev/null +++ b/examples/classification/configs/qa_s06.in @@ -0,0 +1,14 @@ +&config +notrace_passing = 1 +ntestpart = 100000 +num_surf = 1 +sbeg = 0.6 +npoiper2 = 128 +netcdffile = 'wout_henneberg_qa.nc' +isw_field_type = 2 +trace_time = 1.0d0 +tcut = 1d-1 +class_plot = .False. +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qa_volume.in b/examples/classification/configs/qa_volume.in new file mode 100644 index 00000000..16a5e54b --- /dev/null +++ b/examples/classification/configs/qa_volume.in @@ -0,0 +1,15 @@ +&config +notrace_passing = 1 +ntestpart = 100000 +num_surf = 0 +sbeg = 0.2 +npoiper2 = 128 +netcdffile = 'wout_henneberg_qa.nc' +isw_field_type = 2 +trace_time = 1.5d-2 +tcut = 1d-2 +class_plot = .True. +cut_in_per = 0d0 +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qh_fig8.in b/examples/classification/configs/qh_fig8.in new file mode 100644 index 00000000..62e7bdb9 --- /dev/null +++ b/examples/classification/configs/qh_fig8.in @@ -0,0 +1,15 @@ +&config +notrace_passing = 0 +ntestpart = 10000 +num_surf = 1 +sbeg = 0.6 +npoiper2 = 128 +netcdffile = 'wout_qh_8_7.nc' +isw_field_type = 2 +trace_time = 1.0d0 +tcut = 1d-1 +class_plot = .True. +cut_in_per = 0.5d0 +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qh_s03.in b/examples/classification/configs/qh_s03.in new file mode 100644 index 00000000..2b7c797c --- /dev/null +++ b/examples/classification/configs/qh_s03.in @@ -0,0 +1,14 @@ +&config +notrace_passing = 1 +ntestpart = 100000 +num_surf = 1 +sbeg = 0.3 +npoiper2 = 128 +netcdffile = 'wout_qh_8_7.nc' +isw_field_type = 2 +trace_time = 1.0d0 +tcut = 1d-1 +class_plot = .False. +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qh_s06.in b/examples/classification/configs/qh_s06.in new file mode 100644 index 00000000..31442702 --- /dev/null +++ b/examples/classification/configs/qh_s06.in @@ -0,0 +1,14 @@ +&config +notrace_passing = 1 +ntestpart = 100000 +num_surf = 1 +sbeg = 0.6 +npoiper2 = 128 +netcdffile = 'wout_qh_8_7.nc' +isw_field_type = 2 +trace_time = 1.0d0 +tcut = 1d-1 +class_plot = .False. +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qh_volume.in b/examples/classification/configs/qh_volume.in new file mode 100644 index 00000000..537bd5ae --- /dev/null +++ b/examples/classification/configs/qh_volume.in @@ -0,0 +1,15 @@ +&config +notrace_passing = 1 +ntestpart = 100000 +num_surf = 0 +sbeg = 0.2 +npoiper2 = 128 +netcdffile = 'wout_qh_8_7.nc' +isw_field_type = 2 +trace_time = 1.5d-2 +tcut = 1d-2 +class_plot = .True. +cut_in_per = 0d0 +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qi_fig8.in b/examples/classification/configs/qi_fig8.in new file mode 100644 index 00000000..d1600bed --- /dev/null +++ b/examples/classification/configs/qi_fig8.in @@ -0,0 +1,15 @@ +&config +notrace_passing = 0 +ntestpart = 10000 +num_surf = 1 +sbeg = 0.6 +npoiper2 = 128 +netcdffile = 'wout_23_1900_fix_bdry.nc' +isw_field_type = 2 +trace_time = 1.0d0 +tcut = 1d-1 +class_plot = .True. +cut_in_per = 0.5d0 +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qi_s03.in b/examples/classification/configs/qi_s03.in new file mode 100644 index 00000000..712ea192 --- /dev/null +++ b/examples/classification/configs/qi_s03.in @@ -0,0 +1,14 @@ +&config +notrace_passing = 1 +ntestpart = 100000 +num_surf = 1 +sbeg = 0.3 +npoiper2 = 128 +netcdffile = 'wout_23_1900_fix_bdry.nc' +isw_field_type = 2 +trace_time = 1.0d0 +tcut = 1d-1 +class_plot = .False. +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qi_s06.in b/examples/classification/configs/qi_s06.in new file mode 100644 index 00000000..66cb4f94 --- /dev/null +++ b/examples/classification/configs/qi_s06.in @@ -0,0 +1,14 @@ +&config +notrace_passing = 1 +ntestpart = 100000 +num_surf = 1 +sbeg = 0.6 +npoiper2 = 128 +netcdffile = 'wout_23_1900_fix_bdry.nc' +isw_field_type = 2 +trace_time = 1.0d0 +tcut = 1d-1 +class_plot = .False. +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/configs/qi_volume.in b/examples/classification/configs/qi_volume.in new file mode 100644 index 00000000..a004b664 --- /dev/null +++ b/examples/classification/configs/qi_volume.in @@ -0,0 +1,15 @@ +&config +notrace_passing = 1 +ntestpart = 100000 +num_surf = 0 +sbeg = 0.2 +npoiper2 = 128 +netcdffile = 'wout_23_1900_fix_bdry.nc' +isw_field_type = 2 +trace_time = 1.5d-2 +tcut = 1d-2 +class_plot = .True. +cut_in_per = 0d0 +fast_class = .True. +deterministic = .True. +/ diff --git a/examples/classification/plot_classification.py b/examples/classification/plot_classification.py index f0a63ea5..a6efeec6 100644 --- a/examples/classification/plot_classification.py +++ b/examples/classification/plot_classification.py @@ -1,63 +1,98 @@ -import numpy as np -import matplotlib.pyplot as plt -import os - -def doplot_inner(prompt, regular, stochastic, bminmax, outfile, arrows): - tot1=prompt+regular+stochastic - wpr=-1. - wst=0. - wre=1. - cla1=wre*regular+wst*stochastic+wpr*prompt - cla1[tot1==0]=np.nan - tot1[tot1==0]=np.nan - cla1=cla1/tot1 - - bmin=min(bminmax[:,1]) - bmax=max(bminmax[:,2]) - jpmin=bmin/bmax - print("Minimum J_perp: ", jpmin) - - plt.figure(figsize=(5,4)) - plt.imshow(cla1.T, origin='lower', extent=[0,1,0,1], aspect='auto', - vmin=-1.0, vmax=1.0, clim=(-1.0,1.0)) - plt.plot(bminmax[:,0],bmin/bminmax[:,1],'k') - plt.plot(bminmax[:,0],bmin/bminmax[:,2],'k--') - plt.plot([.25, .25],[0.0, bmin/bminmax[250,1]], 'k', lw=0.75) - for a in arrows: - plt.arrow(a[0], a[1], a[2], a[3], length_includes_head=True, - head_width=0.01, head_length=0.04, ec='w', fc='w') - plt.xticks([0,.25,.5,.75,1]) - plt.ylim(jpmin, 1) - plt.xlabel('Normalized toroidal flux $s$') - plt.ylabel('Perpendicular invariant $J_\perp$') - - plt.tight_layout() - -def doplot(basedir, outfile1, outfile2, arrows): - # J_parallel regular/chaotic classification - prompt1 = np.loadtxt(os.path.join(basedir, 'prompt1.dat')) - regular1 = np.loadtxt(os.path.join(basedir, 'regular1.dat')) - stochastic1 = np.loadtxt(os.path.join(basedir, 'stochastic1.dat')) - - # Topological ideal/non-ideal classification - prompt2 = np.loadtxt(os.path.join(basedir, 'prompt2.dat')) - regular2 = np.loadtxt(os.path.join(basedir, 'regular2.dat')) - stochastic2 = np.loadtxt(os.path.join(basedir, 'stochastic2.dat')) - - bminmax = np.loadtxt(os.path.join(basedir, 'bminmax.dat')) - - doplot_inner(prompt1, regular1, stochastic1, bminmax, outfile1, []) - plt.title(r'$J_\parallel$ classifier') - cb = plt.colorbar() - cb.set_ticks([-1, 0, 1]) - cb.set_ticklabels(['Prompt losses', 'Non-ideal', 'Prompt']) - plt.savefig(outfile1) - doplot_inner(prompt2, regular2, stochastic2, bminmax, outfile2, arrows) - plt.title('Topological classifier') - # Set colorbar tick labels - cb = plt.colorbar() - cb.set_ticks([-1, 0, 1]) - cb.set_ticklabels(['Prompt losses', 'Non-ideal', 'Ideal']) - plt.savefig(outfile2) - -doplot('.', 'class_jpar.pdf', 'class_ideal.pdf', []) +import numpy as np +import matplotlib.pyplot as plt + + +def bin_classification(s, perp_inv, iclass, ns=50, nperp=100): + """Bin particles into (s, perp_inv) grid by classification type. + + Returns prompt, regular, stochastic arrays of shape (nperp, ns). + """ + hs = 1.0 / ns + pmin = 0.0 + pmax = np.max(perp_inv) + hp = (pmax - pmin) / nperp + + prompt = np.zeros((nperp, ns)) + regular = np.zeros((nperp, ns)) + stochastic = np.zeros((nperp, ns)) + + for ipart in range(len(s)): + i = min(ns, max(1, int(np.ceil(s[ipart] / hs)))) - 1 + k = min(nperp, max(1, int(np.ceil(perp_inv[ipart] / hp)))) - 1 + + if iclass[ipart] == 0: + prompt[k, i] += 1.0 + elif iclass[ipart] == 1: + regular[k, i] += 1.0 + elif iclass[ipart] == 2: + stochastic[k, i] += 1.0 + + return prompt, regular, stochastic + + +def doplot_inner(prompt, regular, stochastic, bminmax): + tot = prompt + regular + stochastic + wpr = -1.0 + wst = 0.0 + wre = 1.0 + cla = wre * regular + wst * stochastic + wpr * prompt + cla[tot == 0] = np.nan + tot[tot == 0] = np.nan + cla = cla / tot + + bmin_global = np.min(bminmax[:, 1]) + bmax_global = np.max(bminmax[:, 2]) + jpmin = bmin_global / bmax_global + + plt.figure(figsize=(5, 4)) + plt.imshow(cla, origin='lower', extent=[0, 1, 0, 1], aspect='auto', + vmin=-1.0, vmax=1.0, clim=(-1.0, 1.0)) + plt.plot(bminmax[:, 0], bmin_global / bminmax[:, 1], 'k') + plt.plot(bminmax[:, 0], bmin_global / bminmax[:, 2], 'k--') + + # Vertical line at s=0.25 + idx_025 = np.argmin(np.abs(bminmax[:, 0] - 0.25)) + plt.plot([0.25, 0.25], [0.0, bmin_global / bminmax[idx_025, 1]], + 'k', lw=0.75) + + plt.xticks([0, 0.25, 0.5, 0.75, 1]) + plt.ylim(jpmin, 1) + plt.xlabel(r'Normalized toroidal flux $s$') + plt.ylabel(r'Perpendicular invariant $J_\perp$') + plt.tight_layout() + + +def main(): + # Read class_parts.dat: columns are + # ipart, s, perp_inv, iclass_jpar, iclass_ideal, iclass_fractal + data = np.loadtxt('class_parts.dat') + s = data[:, 1] + perp_inv = data[:, 2] + icl_jpar = data[:, 3].astype(int) + icl_ideal = data[:, 4].astype(int) + + bminmax = np.loadtxt('bminmax.dat') + + # J_parallel classification + prompt1, regular1, stochastic1 = bin_classification( + s, perp_inv, icl_jpar) + doplot_inner(prompt1, regular1, stochastic1, bminmax) + plt.title(r'$J_\parallel$ classifier') + cb = plt.colorbar() + cb.set_ticks([-1, 0, 1]) + cb.set_ticklabels(['Prompt losses', 'Non-ideal', 'Regular']) + plt.savefig('class_jpar.pdf') + + # Topological (ideal/non-ideal) classification + prompt2, regular2, stochastic2 = bin_classification( + s, perp_inv, icl_ideal) + doplot_inner(prompt2, regular2, stochastic2, bminmax) + plt.title('Topological classifier') + cb = plt.colorbar() + cb.set_ticks([-1, 0, 1]) + cb.set_ticklabels(['Prompt losses', 'Non-ideal', 'Ideal']) + plt.savefig('class_ideal.pdf') + + +if __name__ == '__main__': + main() diff --git a/examples/classification/plot_paper_results.py b/examples/classification/plot_paper_results.py new file mode 100644 index 00000000..cfe1d610 --- /dev/null +++ b/examples/classification/plot_paper_results.py @@ -0,0 +1,504 @@ +#!/usr/bin/env python3 +""" +Reproduce all classification figures from the JPP 2020 paper: + Albert, Kasilov, Kernbichler, "Accelerated methods for direct computation + of fusion alpha particle losses within stellarator optimization", + J. Plasma Phys. 86, 815860201 (2020). + +Uses the actual paper equilibria: + - QI: Drevlak et al. 2014 (wout_23_1900_fix_bdry.nc, nfp=5) + - QH: Drevlak et al. 2018 (wout_qh_8_7.nc, nfp=5) + - QA: Henneberg et al. 2019 (wout_henneberg_qa.nc, nfp=2) + +Generates: + - Figure 5: QI loss plots at s=0.3 and s=0.6 + - Figure 6: QH loss plots at s=0.3 and s=0.6 + - Figure 7: QA loss plots at s=0.3 and s=0.6 + - Figure 8: Orbit classification in (theta, v_par/v) for all 3 configs + - Table 1: Regular orbit fractions for trapped/passing regions + - Volume classification plots for QH, QI, QA + +Usage: + python plot_paper_results.py [base_dir] + + base_dir defaults to /tmp/simple_classification +""" + +import sys +from pathlib import Path + +import numpy as np + +# Import pysimple.plotting without requiring compiled backend +_plotting_path = ( + Path(__file__).resolve().parent.parent.parent / "python" / "pysimple" / "plotting.py" +) +import importlib.util + +_spec = importlib.util.spec_from_file_location("pysimple.plotting", _plotting_path) +_plotting = importlib.util.module_from_spec(_spec) +sys.modules["pysimple.plotting"] = _plotting +_spec.loader.exec_module(_plotting) + +load_loss_data = _plotting.load_loss_data + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from scipy.stats import gaussian_kde + +# Configuration labels and directory name mappings +CONFIGS = { + "QI": {"label": "QI (Drevlak 2014)", "s03": "qi_s03", "s06": "qi_s06", "fig8": "qi_fig8"}, + "QH": {"label": "QH (Drevlak 2018)", "s03": "qh_s03", "s06": "qh_s06", "fig8": "qh_fig8"}, + "QA": {"label": "QA (Henneberg 2019)", "s03": "qa_s03", "s06": "qa_s06", "fig8": "qa_fig8"}, +} +FIGURE_NUMBERS = {"QI": 5, "QH": 6, "QA": 7} + + +def plot_loss_panel(ax, data, title): + """Plot a single loss panel: KDE density + dots + confined fraction. + + Matches the style of Figures 5-7 in the paper: black dots for lost + particles in (log10(t), trapping parameter) space with blue KDE shading, + a red confined fraction curve on the right axis, and error bands. + """ + lost = data.lost_mask + tlost = np.abs(data.loss_times[lost]) + trap_par = data.trap_parameter[lost] + t_threshold = 0.99 * data.trace_time + kde_mask = tlost < t_threshold + + # KDE density shading + if np.sum(kde_mask) > 10: + log_tlost = np.log10(tlost[kde_mask]) + trap_kde = trap_par[kde_mask] + try: + kde = gaussian_kde([log_tlost, trap_kde]) + xg, yg = np.mgrid[-5:0:100j, -2:1:150j] + positions = np.vstack([xg.ravel(), yg.ravel()]) + density = np.reshape(kde(positions).T, xg.shape) + ax.imshow( + np.rot90(density), cmap=plt.cm.Blues, + extent=[-5, 0, -2, 1], aspect="auto", + ) + except np.linalg.LinAlgError: + pass + + # Lost particle dots + if np.sum(lost) > 0: + ax.plot(np.log10(tlost), trap_par, "k,", markersize=0.5, alpha=0.5) + + ax.axhline(y=0, color="k", linestyle="--", linewidth=0.8) + ax.set_ylim([-2.0, 1.0]) + ax.set_xlim([-5.0, 0.0]) + ax.set_ylabel(r"trapping parameter $\theta_\mathrm{trap}$") + ax.set_xlabel(r"$t$ / s") + ax.set_xticks([-5, -4, -3, -2, -1, 0]) + ax.set_xticklabels( + [r"$10^{-5}$", r"$10^{-4}$", r"$10^{-3}$", + r"$10^{-2}$", r"$10^{-1}$", r"$1$"] + ) + ax.set_title(title) + + # Confined fraction on twin axis + ax2 = ax.twinx() + n = data.n_particles + t_grid = np.logspace(-5, 0, 200) + fc = np.array([np.sum(np.abs(data.loss_times) >= t) / n for t in t_grid]) + sigma = np.sqrt(fc * (1 - fc) / n) + + ax2.plot(np.log10(t_grid), fc, color="tab:red", lw=1.5) + ax2.fill_between( + np.log10(t_grid), fc - 1.96 * sigma, fc + 1.96 * sigma, + color="tab:red", alpha=0.15, + ) + ax2.set_ylim([0.6, 1.05]) + ax2.set_ylabel(r"confined fraction $f_c$", color="tab:red") + ax2.tick_params("y", colors="tab:red") + + +def plot_loss_figure(data_s03, data_s06, config_key, output_path): + """Generate a two-panel loss figure (Figures 5-7 style).""" + fig_num = FIGURE_NUMBERS[config_key] + label = CONFIGS[config_key]["label"] + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + + plot_loss_panel(axes[0], data_s03, f"({chr(97)}) $s = 0.3$") + plot_loss_panel(axes[1], data_s06, f"({chr(98)}) $s = 0.6$") + + fig.suptitle(f"Figure {fig_num}: {label}", fontsize=13, y=1.02) + fig.tight_layout() + fig.savefig(str(output_path), dpi=150, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {output_path}") + + +def _read_namelist_param(simple_in, param_name, default): + """Read a single parameter value from a simple.in namelist file.""" + if not simple_in.exists(): + return default + for line in simple_in.read_text().splitlines(): + if param_name in line and "=" in line: + if param_name == "tcut" and "ntcut" in line: + continue + val = line.split("=")[1].split("!")[0].strip().rstrip(",") + return float(val.replace("d", "e")) + return default + + +def _classify_particles(fig8_dir): + """Load and classify particles from a Fig 8 run directory. + + Returns (theta, pitch, category) where category is: + 0 = regular (confined), 1 = early loss, 2 = chaotic/late loss + Returns None if data files are missing. + """ + for required in ("class_parts.dat", "times_lost.dat", "start.dat"): + if not (fig8_dir / required).exists(): + return None + + class_data = np.loadtxt(fig8_dir / "class_parts.dat") + iclass_jpar = class_data[:, 3].astype(int) + + tl_data = np.loadtxt(fig8_dir / "times_lost.dat") + loss_times = tl_data[:, 1] + + start_data = np.loadtxt(fig8_dir / "start.dat") + theta = np.mod(start_data[:, 1], 2 * np.pi) / np.pi + pitch = start_data[:, 4] + + simple_in = fig8_dir / "simple.in" + trace_time = _read_namelist_param(simple_in, "trace_time", 1.0) + tcut = _read_namelist_param(simple_in, "tcut", 0.1) + + regular = iclass_jpar == 1 + prompt_loss = iclass_jpar == 0 + chaotic = iclass_jpar == 2 + + lost_early = (loss_times > 0) & (loss_times < tcut) + lost_late = (loss_times > 0) & (loss_times >= tcut) & (loss_times < trace_time) + confined = (loss_times < 0) | (loss_times >= trace_time) + + # Category: 0=regular, 1=early loss, 2=chaotic + category = np.full(len(theta), -1, dtype=int) + category[regular] = 0 + category[prompt_loss | (chaotic & lost_early)] = 1 + category[chaotic & (lost_late | confined)] = 2 + # Regular but lost -- still mark as regular for background + category[regular & ~confined] = 0 + + return theta, pitch, category + + +def _plot_fig8_panel(ax, fig8_dir, panel_label): + """Plot a single Fig 8 panel with dense imshow background and overlay markers. + + Uses a coarse grid with Gaussian smoothing so the background fills the + entire (theta/pi, v_par/v) plane -- the "brazilian flag" look from the + paper. Colors: green=regular, blue=early loss, yellow=chaotic. + Overlay markers: circles for early losses, crosses for chaotic. + """ + from scipy.ndimage import gaussian_filter + + result = _classify_particles(fig8_dir) + if result is None: + ax.set_title(f"{panel_label} -- no data") + return + + theta, pitch, category = result + + n_theta, n_pitch = 50, 50 + theta_edges = np.linspace(0, 2, n_theta + 1) + pitch_edges = np.linspace(-1, 1, n_pitch + 1) + + # Count particles of each type in each grid cell + counts = np.zeros((3, n_pitch, n_theta)) + for cat_id in range(3): + mask = category == cat_id + if np.any(mask): + h, _, _ = np.histogram2d( + pitch[mask], theta[mask], + bins=[pitch_edges, theta_edges], + ) + counts[cat_id] = h + + # Smooth counts to fill gaps and produce continuous colored regions + sigma = 1.5 + smoothed = np.array([gaussian_filter(counts[i], sigma) for i in range(3)]) + total_smooth = smoothed.sum(axis=0) + + # Build RGB image from smoothed fractional composition + color_map = np.array([ + [0.2, 0.7, 0.2], # green (regular) + [0.2, 0.4, 0.9], # blue (early loss) + [0.9, 0.8, 0.1], # yellow (chaotic) + ]) + rgb = np.ones((n_pitch, n_theta, 3)) + filled = total_smooth > 1e-10 + for cat_id in range(3): + frac = np.where(filled, smoothed[cat_id] / np.where(filled, total_smooth, 1.0), 0.0) + for ch in range(3): + rgb[:, :, ch] = np.where( + filled, + rgb[:, :, ch] - frac * (1.0 - color_map[cat_id, ch]), + rgb[:, :, ch], + ) + rgb = np.clip(rgb, 0, 1) + + ax.imshow( + rgb, origin="lower", extent=[0, 2, -1, 1], aspect="auto", + interpolation="bilinear", + ) + + # Overlay markers for early losses (circles) and chaotic (crosses) + early = category == 1 + if np.any(early): + ax.scatter( + theta[early], pitch[early], + c="blue", s=2, marker="o", alpha=0.3, edgecolors="none", + label=r"Early loss ($t < t_\mathrm{cut}$)", + ) + + chaotic_mask = category == 2 + if np.any(chaotic_mask): + ax.scatter( + theta[chaotic_mask], pitch[chaotic_mask], + c="red", s=4, marker="x", alpha=0.4, linewidths=0.5, + label="Chaotic", + ) + + # Trapped-passing boundary + ax.axhline(y=0, color="white", linewidth=1.5, alpha=0.8) + + ax.set_xlabel(r"$\vartheta / \pi$") + ax.set_xlim([0, 2]) + ax.set_ylim([-1, 1]) + ax.set_title(panel_label) + + +def plot_fig8(base_dir, output_path): + """Generate Figure 8: orbit classification in (theta, v_par/v) space. + + Three panels (a, b, c) for QI, QH, QA at s=0.6, phi = phi_n/2. + Uses a dense colored grid (imshow) for the background with overlay markers. + """ + fig, axes = plt.subplots(1, 3, figsize=(14, 4.5)) + panel_labels = ["(a) QI", "(b) QH", "(c) QA"] + + for ax, config_key, panel_label in zip(axes, ["QI", "QH", "QA"], panel_labels): + fig8_dir = base_dir / CONFIGS[config_key]["fig8"] + _plot_fig8_panel(ax, fig8_dir, panel_label) + + # Only the leftmost panel gets a y-label + axes[0].set_ylabel(r"$v_\parallel / v$") + + # Shared legend from whichever panel has handles + handles, labels = [], [] + for ax in axes: + h, l = ax.get_legend_handles_labels() + for hi, li in zip(h, l): + if li not in labels: + handles.append(hi) + labels.append(li) + if handles: + fig.legend(handles, labels, loc="lower center", ncol=4, fontsize=8, + bbox_to_anchor=(0.5, -0.05)) + + fig.suptitle("Figure 8: Orbit classification at $s = 0.6$", fontsize=13, y=1.02) + fig.tight_layout() + fig.savefig(str(output_path), dpi=150, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {output_path}") + + +def plot_volume_classification(data_dir, output_path, config_label=""): + """Volume classification plot: (s, J_perp) colored by orbit type. + + Bins trapped particles only (iclass != 1 or trap_par >= 0) into a + (s, J_perp) grid. Passing particles are shown as white background. + The trapped-passing and deeply-trapped boundaries from bminmax.dat + are overlaid. + """ + data_dir = Path(data_dir) + + class_data = np.loadtxt(data_dir / "class_parts.dat") + s = class_data[:, 1] + perp_inv = class_data[:, 2] + icl_jpar = class_data[:, 3].astype(int) + + tl_data = np.loadtxt(data_dir / "times_lost.dat") + trap_par = tl_data[:, 2] + + bminmax = np.loadtxt(data_dir / "bminmax.dat") + + ns, nperp = 50, 100 + hs = 1.0 / ns + pmax = np.max(perp_inv) + hp = pmax / nperp + + prompt = np.zeros((nperp, ns)) + regular = np.zeros((nperp, ns)) + stochastic = np.zeros((nperp, ns)) + + for ipart in range(len(s)): + i = min(ns, max(1, int(np.ceil(s[ipart] / hs)))) - 1 + k = min(nperp, max(1, int(np.ceil(perp_inv[ipart] / hp)))) - 1 + is_trapped = trap_par[ipart] >= 0 + if not is_trapped: + continue + if icl_jpar[ipart] == 0: + prompt[k, i] += 1.0 + elif icl_jpar[ipart] == 1: + regular[k, i] += 1.0 + elif icl_jpar[ipart] == 2: + stochastic[k, i] += 1.0 + + tot = prompt + regular + stochastic + cla = regular - prompt + cla[tot == 0] = np.nan + tot[tot == 0] = np.nan + cla = cla / tot + + bmin_global = np.min(bminmax[:, 1]) + + fig, ax = plt.subplots(figsize=(6, 5)) + im = ax.imshow( + cla, origin="lower", extent=[0, 1, 0, 1], aspect="auto", + vmin=-1.0, vmax=1.0, + ) + ax.plot(bminmax[:, 0], bmin_global / bminmax[:, 1], "k", lw=1.5) + ax.plot(bminmax[:, 0], bmin_global / bminmax[:, 2], "k--", lw=1.0) + + jpmin = bmin_global / np.max(bminmax[:, 2]) + ax.set_ylim(jpmin, 1) + ax.set_xlabel(r"Normalized toroidal flux $s$") + ax.set_ylabel(r"Perpendicular invariant $J_\perp$") + title_suffix = f" -- {config_label}" if config_label else "" + ax.set_title(f"Volume classification{title_suffix}") + + cb = plt.colorbar(im, ax=ax) + cb.set_ticks([-1, 0, 1]) + cb.set_ticklabels(["Prompt losses", "Chaotic", "Regular"]) + + fig.tight_layout() + fig.savefig(str(output_path), dpi=150, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {output_path}") + + +def print_table1(base_dir): + """Print Table 1: fractions of regular orbits in trapped/passing regions.""" + print("\n Table 1: Fractions of regular orbits") + print(" " + "-" * 60) + print(f" {'Type':>6s} {'Surface s':>10s} {'Regular trapped':>16s} {'Regular passing':>16s}") + print(" " + "-" * 60) + + for config_key in ["QI", "QH", "QA"]: + for surface, dir_name in [("0.3", CONFIGS[config_key]["s03"]), + ("0.6", CONFIGS[config_key]["s06"])]: + run_dir = base_dir / dir_name + class_file = run_dir / "class_parts.dat" + tl_file = run_dir / "times_lost.dat" + + if not class_file.exists() or not tl_file.exists(): + print(f" {config_key:>6s} {surface:>10s} {'N/A':>16s} {'N/A':>16s}") + continue + + class_data = np.loadtxt(class_file) + icl_jpar = class_data[:, 3].astype(int) + + tl_data = np.loadtxt(tl_file) + trap_par = tl_data[:, 2] + + trapped = trap_par >= 0 + passing = trap_par < 0 + + n_trapped = np.sum(trapped) + n_passing = np.sum(passing) + + regular_trapped = np.sum((icl_jpar == 1) & trapped) + regular_passing = np.sum((icl_jpar == 1) & passing) + + frac_trapped = regular_trapped / n_trapped if n_trapped > 0 else 0.0 + frac_passing = regular_passing / n_passing if n_passing > 0 else 0.0 + + print(f" {config_key:>6s} {surface:>10s} {frac_trapped:>16.4f} {frac_passing:>16.4f}") + + print(" " + "-" * 60) + + +def print_loss_summary(data, label): + """Print loss summary for a single-surface run.""" + n = data.n_particles + lost = data.lost_mask + n_lost = np.sum(lost) + fc = 1.0 - n_lost / n + sigma = np.sqrt(fc * (1 - fc) / n) + print(f" {label}: {n} particles, {n_lost} lost, " + f"f_c = {fc:.4f} +/- {1.96 * sigma:.4f}") + + +def main(): + if len(sys.argv) > 1: + base_dir = Path(sys.argv[1]) + else: + base_dir = Path("/tmp/simple_classification") + + output_dir = base_dir / "plots" + output_dir.mkdir(exist_ok=True) + + # Load data for all configurations + datasets = {} + for config_key in ["QI", "QH", "QA"]: + for surface in ["s03", "s06"]: + dir_name = CONFIGS[config_key][surface] + run_dir = base_dir / dir_name + tl_file = run_dir / "times_lost.dat" + key = f"{config_key}_{surface}" + if tl_file.exists(): + datasets[key] = load_loss_data(run_dir) + s_label = "0.3" if surface == "s03" else "0.6" + print_loss_summary(datasets[key], f"{config_key} s={s_label}") + + # Figures 5-7: loss plots for each configuration + for config_key, fig_num in FIGURE_NUMBERS.items(): + key_s03 = f"{config_key}_s03" + key_s06 = f"{config_key}_s06" + if key_s03 in datasets and key_s06 in datasets: + print(f"\nGenerating Figure {fig_num} ({config_key} loss plots)...") + plot_loss_figure( + datasets[key_s03], datasets[key_s06], config_key, + output_dir / f"fig{fig_num}_losses_{config_key.lower()}.pdf", + ) + + # Figure 8: orbit classification + have_any_fig8 = any( + (base_dir / CONFIGS[ck]["fig8"] / "class_parts.dat").exists() + for ck in ["QI", "QH", "QA"] + ) + if have_any_fig8: + print("\nGenerating Figure 8 (orbit classification)...") + plot_fig8(base_dir, output_dir / "fig8_classification.pdf") + + # Volume classification for all 3 configs + volume_configs = [ + ("QH", "qh_volume"), + ("QI", "qi_volume"), + ("QA", "qa_volume"), + ] + for vol_label, vol_name in volume_configs: + vol_dir = base_dir / vol_name + if (vol_dir / "class_parts.dat").exists(): + out_name = f"volume_classification_{vol_label.lower()}.pdf" + print(f"\nGenerating volume classification plot ({vol_label})...") + plot_volume_classification(vol_dir, output_dir / out_name, vol_label) + + # Table 1 + print_table1(base_dir) + + print(f"\nAll plots saved to: {output_dir}") + + +if __name__ == "__main__": + main() diff --git a/examples/classification/postprocess_class.f90 b/examples/classification/postprocess_class.f90 deleted file mode 100644 index 70beee0f..00000000 --- a/examples/classification/postprocess_class.f90 +++ /dev/null @@ -1,100 +0,0 @@ -! - implicit none - ! - integer, parameter :: ns=50,nperp=100 - integer :: npart,i,ipart,j,k - double precision :: dummy,pmax,pmin,hs,hp - double precision, dimension(nperp,ns,2) :: promp,regul,stoch - integer, dimension(:), allocatable :: icl_jpa,icl_rec - double precision, dimension(:), allocatable :: s,perp_inv - ! - npart=0 - open(1,file='class_parts.dat') - do - read (1,*,end=1) j - npart=npart+1 - enddo - 1 close(1) - print *,'npart = ',npart - ! - allocate(icl_jpa(npart),icl_rec(npart),s(npart),perp_inv(npart)) - ! - open(1,file='class_parts.dat') - do ipart=1,npart - read (1,*) j,s(ipart),perp_inv(ipart),icl_jpa(ipart),icl_rec(ipart) - enddo - close(1) - ! - hs=1.d0/dble(ns) - ! - pmin=0.d0 - pmax=maxval(perp_inv) - promp=0.d0 - regul=0.d0 - stoch=0.d0 - ! - hp=(pmax-pmin)/dble(nperp) - ! - do ipart=1,npart - i=ceiling(s(ipart)/hs) - i=min(ns,max(1,i)) - ! - k=ceiling(perp_inv(ipart)/hp) - k=min(nperp,max(1,k)) - ! - select case(icl_jpa(ipart)) - case(0) - promp(k,i,1)=promp(k,i,1)+1.d0 - case(1) - regul(k,i,1)=regul(k,i,1)+1.d0 - case(2) - stoch(k,i,1)=stoch(k,i,1)+1.d0 - end select - ! - select case(icl_rec(ipart)) - case(0) - promp(k,i,2)=promp(k,i,2)+1.d0 - case(1) - regul(k,i,2)=regul(k,i,2)+1.d0 - case(2) - stoch(k,i,2)=stoch(k,i,2)+1.d0 - end select - enddo - ! - open(1,file='prompt1.dat') - do i=1,ns - write(1,*) promp(:,i,1) - enddo - close(1) - ! - open(1,file='prompt2.dat') - do i=1,ns - write(1,*) promp(:,i,2) - enddo - close(1) - ! - open(1,file='regular1.dat') - do i=1,ns - write(1,*) regul(:,i,1) - enddo - close(1) - ! - open(1,file='regular2.dat') - do i=1,ns - write(1,*) regul(:,i,2) - enddo - close(1) - ! - open(1,file='stochastic1.dat') - do i=1,ns - write(1,*) stoch(:,i,1) - enddo - close(1) - ! - open(1,file='stochastic2.dat') - do i=1,ns - write(1,*) stoch(:,i,2) - enddo - close(1) - ! - end diff --git a/examples/classification/run_all.sh b/examples/classification/run_all.sh new file mode 100755 index 00000000..12e9d1f8 --- /dev/null +++ b/examples/classification/run_all.sh @@ -0,0 +1,62 @@ +#!/bin/bash +# Run all classification cases from the JPP 2020 paper. +# Equilibrium files are downloaded automatically on first run. +set -euo pipefail + +SIMPLE_X="$(cd "$(dirname "$0")/../.." && pwd)/build/simple.x" +EXAMPLE_DIR="$(cd "$(dirname "$0")" && pwd)" +BASE_DIR="${1:-/tmp/simple_classification}" + +if [ ! -x "$SIMPLE_X" ]; then + echo "ERROR: simple.x not found. Build with: make" + exit 1 +fi + +mkdir -p "$BASE_DIR" + +# Download equilibrium files if not already present +download_wout() { + local name="$1" url="$2" dest="$BASE_DIR/$name" + if [ ! -f "$dest" ]; then + echo "Downloading $name ..." + wget -q -O "$dest" "$url" + fi +} + +download_wout wout_qh_8_7.nc \ + "https://github.com/itpplasma/proxima-simple-classification/raw/main/data/wout_qh_8_7.nc" +download_wout wout_23_1900_fix_bdry.nc \ + "https://github.com/itpplasma/proxima-simple-classification/raw/main/data/wout_23_1900_fix_bdry.nc" +download_wout wout_henneberg_qa.nc \ + "https://github.com/itpplasma/proxima-simple-classification/raw/main/data/wout_henneberg_qa.nc" + +run_case() { + local case_name="$1" + local dir="$BASE_DIR/$case_name" + + mkdir -p "$dir" + cp "$EXAMPLE_DIR/configs/${case_name}.in" "$dir/simple.in" + + for wout in "$BASE_DIR"/wout_*.nc; do + ln -sf "$wout" "$dir/" + done + + echo "=== Running $case_name ===" + local t0; t0=$(date +%s) + (cd "$dir" && "$SIMPLE_X") + echo "=== $case_name done in $(($(date +%s) - t0))s ===" +} + +if [ $# -gt 1 ]; then + shift # skip base_dir + for case_name in "$@"; do run_case "$case_name"; done +else + for case_name in \ + qi_s03 qi_s06 qh_s03 qh_s06 qa_s03 qa_s06 \ + qi_fig8 qh_fig8 qa_fig8 \ + qi_volume qh_volume qa_volume; do + run_case "$case_name" + done +fi + +echo "All runs complete. Results in $BASE_DIR" diff --git a/examples/classification/simple.in b/examples/classification/simple.in index b6124283..fad5dd1a 100644 --- a/examples/classification/simple.in +++ b/examples/classification/simple.in @@ -1,15 +1,15 @@ &config -notrace_passing = 1 ! skip tracing passing prts if notrace_passing=1 -ntestpart = 5000 ! number of test particles -num_surf = 0 ! number of start surfaces, start points everywhere if 0 -sbeg = 0.2 ! starting surface for some initializations -npoiper2 = 128 ! integration steps per field period -netcdffile = 'wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc' ! name of VMEC file in NETCDF format -isw_field_type = 2 ! -1: Testing, 0: Canonical, 1: VMEC, 2: Boozer -trace_time = 1.5d-2 ! tracing time -tcut = 1d-2 ! time when to do cut for classification, usually 1d-1, or -1 if no cuts desired -class_plot = .True. ! write starting points at phi=const cut for classification plot (.True./.False.) -cut_in_per = 0d0 ! normalized phi-cut position within field period, [0:1], used if class_plot=.True. -fast_class = .True. ! if .True. quit immeadiately after fast classification and don't trace orbits to the end +notrace_passing = 1 +ntestpart = 5000 +num_surf = 0 +sbeg = 0.2 +npoiper2 = 128 +netcdffile = 'wout_qh_8_7.nc' +isw_field_type = 2 +trace_time = 1.5d-2 +tcut = 1d-2 +class_plot = .True. +cut_in_per = 0d0 +fast_class = .True. deterministic = .True. / diff --git a/src/classification.f90 b/src/classification.f90 index 2836bb14..320f21fa 100644 --- a/src/classification.f90 +++ b/src/classification.f90 @@ -146,7 +146,7 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) call magfie(z(1:3),bmod,sqrtg,bder,hcovar,hctrvr,hcurl) !$omp critical - if(num_surf > 1) then + if(num_surf /= 1) then call get_bminmax(z(1),bmin,bmax) endif passing = z(5)**2.gt.1.d0-bmod/bmax @@ -325,7 +325,13 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) ! Store in classification result class_result%jpar = ijpar class_result%topology = ideal - if(fast_class) ierr=ierr_cot + ! fast_class without class_plot: stop tracing regular + ! orbits early (they are confined). With class_plot, + ! continue to ntcut so classification output is written. + if(fast_class .and. .not. class_plot & + .and. ierr_cot /= 0 .and. ijpar == 1) then + regular = .True. + endif ! ! End classification by J_parallel and ideal orbit conditions endif diff --git a/src/find_bminmax.f90 b/src/find_bminmax.f90 index 545a5ee3..4ecd01c7 100644 --- a/src/find_bminmax.f90 +++ b/src/find_bminmax.f90 @@ -155,33 +155,36 @@ end subroutine find_bminmax ! !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! - subroutine get_bminmax(s,bmin,bmax) - ! - use bminmax_mod, only : prop,nsbmnx,hsbmnx,bmin_arr,bmax_arr - ! - implicit none - ! - integer :: k - real(dp) :: s,bmin,bmax,ws,s0 - ! - if(prop) then - prop=.false. - hsbmnx=1.d0/dble(nsbmnx) - do k=0,nsbmnx - s0=max(1.d-8,hsbmnx*dble(k)) - ! - call find_bminmax(s0,bmin_arr(k),bmax_arr(k)) - ! - enddo - endif - ! - ws=s/hsbmnx - k=min(nsbmnx-1,max(0,int(ws))) - ws=ws-dble(k) - ! - bmin=bmin_arr(k)*(1.d0-ws)+bmin_arr(k+1)*ws - bmax=bmax_arr(k)*(1.d0-ws)+bmax_arr(k+1)*ws - ! + subroutine init_bminmax_arrays + use bminmax_mod, only: prop, nsbmnx, hsbmnx, bmin_arr, bmax_arr + + integer :: k + real(dp) :: s0 + + if (.not. prop) return + prop = .false. + hsbmnx = 1.d0/dble(nsbmnx) + do k = 0, nsbmnx + s0 = max(1.d-8, hsbmnx*dble(k)) + call find_bminmax(s0, bmin_arr(k), bmax_arr(k)) + end do + end subroutine init_bminmax_arrays + + subroutine get_bminmax(s, bmin, bmax) + use bminmax_mod, only: prop, nsbmnx, hsbmnx, bmin_arr, bmax_arr + + real(dp), intent(in) :: s + real(dp), intent(out) :: bmin, bmax + integer :: k + real(dp) :: ws + + if (prop) call init_bminmax_arrays + + ws = s/hsbmnx + k = min(nsbmnx - 1, max(0, int(ws))) + ws = ws - dble(k) + bmin = bmin_arr(k)*(1.d0 - ws) + bmin_arr(k + 1)*ws + bmax = bmax_arr(k)*(1.d0 - ws) + bmax_arr(k + 1)*ws end subroutine get_bminmax end module find_bminmax_sub diff --git a/src/params.f90 b/src/params.f90 index 44288722..19887dfb 100644 --- a/src/params.f90 +++ b/src/params.f90 @@ -137,9 +137,6 @@ subroutine read_config(config_file) error stop 'Collisions are incompatible with classification' end if - if (fast_class .and. tcut > 0.0d0) then - error stop 'fast_class and positive tcut are mutually exclusive' - end if end subroutine read_config subroutine params_init diff --git a/src/simple_main.f90 b/src/simple_main.f90 index cd7d8842..d9dda8a7 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -98,6 +98,9 @@ subroutine main call init_starting_surf call print_phase_time('Starting surface initialization completed') + call init_bminmax + call print_phase_time('Bmin/Bmax initialization completed') + call sample_particles call print_phase_time('Particle sampling completed') @@ -426,7 +429,9 @@ subroutine sample_particles use samplers, only: sample, START_FILE if (1 == startmode) then - if ((0d0 < grid_density) .and. (1d0 > grid_density)) then + if (0 == num_surf) then + call sample(zstart, 0.0d0, 1.0d0) + else if ((0d0 < grid_density) .and. (1d0 > grid_density)) then call sample(zstart, grid_density) else call sample(zstart) @@ -514,6 +519,16 @@ subroutine save_starting_points_test(zs) end subroutine save_starting_points_test end subroutine sample_particles_test_field + subroutine init_bminmax + use find_bminmax_sub, only: init_bminmax_arrays + + ! Populate bminmax arrays while magfie is still in VMEC mode. + ! find_bminmax scans (theta, phi) to find field extrema, so it + ! must run in the coordinate system where theta/phi are geometric + ! angles (VMEC), not canonical angles (Boozer/canflux). + call init_bminmax_arrays + end subroutine init_bminmax + subroutine init_counters icounter = 0 ! evaluation counter kpart = 0 @@ -768,7 +783,7 @@ subroutine compute_pitch_angle_params(z, passing, trap_par_, perp_inv_) !$omp critical bmod = compute_bmod(z(1:3)) - if (num_surf > 1) then + if (num_surf /= 1) then call get_bminmax(z(1), bmin, bmax) end if passing = z(5)**2 .gt. 1.d0 - bmod/bmax @@ -862,6 +877,16 @@ subroutine write_output write (1, *) i, zstart(1, i), perp_inv(i), iclass(:, i) end do close (1) + + block + use bminmax_mod, only: nsbmnx, hsbmnx, bmin_arr, bmax_arr + + open (1, file='bminmax.dat', recl=1024) + do i = 0, nsbmnx + write (1, *) hsbmnx*dble(i), bmin_arr(i), bmax_arr(i) + end do + close (1) + end block end if if (output_results_netcdf) then diff --git a/test/golden_record/classifier_combined/simple.in b/test/golden_record/classifier_combined/simple.in new file mode 100644 index 00000000..6ff1800e --- /dev/null +++ b/test/golden_record/classifier_combined/simple.in @@ -0,0 +1,15 @@ +&config +notrace_passing = 1 ! skip tracing passing prts if notrace_passing=1 +ntestpart = 32 ! number of test particles +num_surf = 0 ! number of start surfaces, start points everywhere if 0 +sbeg = 0.2 ! starting surface for some initializations +npoiper2 = 128 ! integration steps per field period +netcdffile = 'wout.nc' ! name of VMEC file in NETCDF format +isw_field_type = 2 ! -1: Testing, 0: Canonical, 1: VMEC, 2: Boozer +trace_time = 1.5d-2 ! tracing time +tcut = 1d-2 ! time when to do cut for classification (fractal dimension) +class_plot = .True. ! write starting points at phi=const cut for classification plot (.True./.False.) +cut_in_per = 0d0 ! normalized phi-cut position within field period, [0:1], used if class_plot=.True. +fast_class = .True. ! if .True. quit immediately after fast classification +deterministic = .True. +/ diff --git a/test/golden_record/compare_golden_results.sh b/test/golden_record/compare_golden_results.sh index 18add9db..f69ddbbb 100755 --- a/test/golden_record/compare_golden_results.sh +++ b/test/golden_record/compare_golden_results.sh @@ -123,8 +123,8 @@ compare_cases() { echo " ✗ FAILED" failed_cases=$((failed_cases + 1)) fi - # Check if this is the classifier_fast case with multiple files - elif [ "$CASE" = "classifier_fast" ]; then + # Check if this is a classifier case with multiple files + elif [ "$CASE" = "classifier_fast" ] || [ "$CASE" = "classifier_combined" ]; then # List of files to compare for classifier_fast (excluding simple.in and wout.nc) # Note: fort.* files are excluded due to non-deterministic ordering in parallel execution CLASSIFIER_FILES="avg_inverse_t_lost.dat class_parts.dat confined_fraction.dat healaxis.dat start.dat times_lost.dat" diff --git a/test/tests/magfie/plot_orbit_chartmap_comparison.py b/test/tests/magfie/plot_orbit_chartmap_comparison.py index a5f58973..acf7357c 100644 --- a/test/tests/magfie/plot_orbit_chartmap_comparison.py +++ b/test/tests/magfie/plot_orbit_chartmap_comparison.py @@ -6,14 +6,19 @@ comparing all four integration paths in Cartesian coordinates. """ -import numpy as np -import matplotlib -matplotlib.use("Agg", force=True) -import matplotlib.pyplot as plt -import netCDF4 as nc import sys from pathlib import Path +try: + import numpy as np + import matplotlib + matplotlib.use("Agg", force=True) + import matplotlib.pyplot as plt + import netCDF4 as nc +except ImportError as e: + print(f"Skipping plot test: {e}") + sys.exit(0) + def load_data(filename): """Load trajectory data from NetCDF file."""