From bd6e1521b06b1620ccac6ac8ac8c9d643fcbfb4a Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Fri, 18 Oct 2024 15:56:17 +0200 Subject: [PATCH] tests: check EC profiles --- scripts/gray_visual.py | 2 +- src/gray_tables.f90 | 2 +- tests/01-ITER/outputs/ec-profiles.48.txt | 2 +- .../outputs/ec-profiles.48.txt | 2 +- tests/03-TCV/outputs/ec-profiles.48.txt | 2 +- tests/04-JT60SA/outputs/ec-profiles.48.txt | 2 +- .../outputs/ec-profiles.48.txt | 2 +- .../outputs/mixed/ec-profiles.48.txt | 2 +- .../outputs/x-mode/ec-profiles.48.txt | 2 +- tests/07-DTT/outputs/ec-profiles.48.txt | 2 +- tests/08-STEP/outputs/ec-profiles.48.txt | 2 +- .../09-DTT-rampup/outputs/ec-profiles.48.txt | 2 +- tests/12-sum-profiles/__init__.py | 40 +++++++++++-- tests/__init__.py | 60 +++++++++++++++++++ 14 files changed, 107 insertions(+), 17 deletions(-) diff --git a/scripts/gray_visual.py b/scripts/gray_visual.py index 3145762..6fbd432 100644 --- a/scripts/gray_visual.py +++ b/scripts/gray_visual.py @@ -351,7 +351,7 @@ def plot_profiles(outputs: Path, axes: [plt.Axes]): slice = profiles[ (profiles['index_rt'] == i) & (profiles['dPdV'] > 1e-15)] - cur_dens.plot(slice['ρ_p'], slice['J_cdb'], **style) + cur_dens.plot(slice['ρ_p'], slice['J_cd'], **style) pow_dens.plot(slice['ρ_p'], slice['dPdV'], **style) slice = profiles[ (profiles['index_rt'] == i) diff --git a/src/gray_tables.f90 b/src/gray_tables.f90 index 52eac3e..e49f4c6 100644 --- a/src/gray_tables.f90 +++ b/src/gray_tables.f90 @@ -47,7 +47,7 @@ contains call tbls%ec_profiles%init( & title='ec-profiles', id=48, active=is_active(params, 48), & labels=[character(64) :: 'ρ_p', 'ρ_t', 'J_φ', & - 'J_cdb', 'dPdV', 'I_cd_inside', & + 'J_cd', 'dPdV', 'I_cd_inside', & 'P_inside', 'index_rt']) call tbls%summary%init( & diff --git a/tests/01-ITER/outputs/ec-profiles.48.txt b/tests/01-ITER/outputs/ec-profiles.48.txt index 179f552..e9d185e 100644 --- a/tests/01-ITER/outputs/ec-profiles.48.txt +++ b/tests/01-ITER/outputs/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 11 16 1800 1.000 0.2000 # COD iwarm ilarm imx ieccd: 2 5 -20 3 # COD ipec nrho istpr istpl: 1 501 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 2.00000000E-003 1.52488048E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 4.00000000E-003 3.04976097E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 diff --git a/tests/02-ITER-half-field/outputs/ec-profiles.48.txt b/tests/02-ITER-half-field/outputs/ec-profiles.48.txt index 35a44ed..ba1f134 100644 --- a/tests/02-ITER-half-field/outputs/ec-profiles.48.txt +++ b/tests/02-ITER-half-field/outputs/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 11 16 2800 1.000 0.2000 # COD iwarm ilarm imx ieccd: 2 5 -20 3 # COD ipec nrho istpr istpl: 1 251 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 4.00000000E-003 3.15628723E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 8.00000000E-003 6.31257445E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 diff --git a/tests/03-TCV/outputs/ec-profiles.48.txt b/tests/03-TCV/outputs/ec-profiles.48.txt index c2c332c..8180451 100644 --- a/tests/03-TCV/outputs/ec-profiles.48.txt +++ b/tests/03-TCV/outputs/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 16 24 2400 1.000 0.5000E-1 # COD iwarm ilarm imx ieccd: 2 5 -20 3 # COD ipec nrho istpr istpl: 1 1001 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 1.00000000E-003 7.42366033E-004 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 2.00000000E-003 1.48473207E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 diff --git a/tests/04-JT60SA/outputs/ec-profiles.48.txt b/tests/04-JT60SA/outputs/ec-profiles.48.txt index f77708d..db71182 100644 --- a/tests/04-JT60SA/outputs/ec-profiles.48.txt +++ b/tests/04-JT60SA/outputs/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 1 1 15000 1.200 0.1000 # COD iwarm ilarm imx ieccd: 2 4 -20 3 # COD ipec nrho istpr istpl: 1 2001 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 5.00000000E-004 4.10134311E-004 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 1.00000000E-003 8.20268623E-004 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 diff --git a/tests/05-JT60SA-startup/outputs/ec-profiles.48.txt b/tests/05-JT60SA-startup/outputs/ec-profiles.48.txt index 1c4154d..fd7ee3c 100644 --- a/tests/05-JT60SA-startup/outputs/ec-profiles.48.txt +++ b/tests/05-JT60SA-startup/outputs/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 1 1 10000 1.000 0.1000 # COD iwarm ilarm imx ieccd: 2 4 -20 3 # COD ipec nrho istpr istpl: 1 2001 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 5.00000000E-004 3.70690524E-004 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 1.00000000E-003 7.41381049E-004 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 diff --git a/tests/06-ITER-startup/outputs/mixed/ec-profiles.48.txt b/tests/06-ITER-startup/outputs/mixed/ec-profiles.48.txt index da4e593..c9eeea6 100644 --- a/tests/06-ITER-startup/outputs/mixed/ec-profiles.48.txt +++ b/tests/06-ITER-startup/outputs/mixed/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 11 16 8000 1.000 0.1000 # COD iwarm ilarm imx ieccd: 2 4 -20 3 # COD ipec nrho istpr istpl: 1 501 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 2.00000000E-003 2.22948719E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 4.00000000E-003 4.45897439E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 diff --git a/tests/06-ITER-startup/outputs/x-mode/ec-profiles.48.txt b/tests/06-ITER-startup/outputs/x-mode/ec-profiles.48.txt index 6ca9064..dbf71af 100644 --- a/tests/06-ITER-startup/outputs/x-mode/ec-profiles.48.txt +++ b/tests/06-ITER-startup/outputs/x-mode/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 11 16 8000 1.000 0.1000 # COD iwarm ilarm imx ieccd: 2 4 -20 3 # COD ipec nrho istpr istpl: 1 501 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 2.00000000E-003 2.22948719E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 4.00000000E-003 4.45897439E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 diff --git a/tests/07-DTT/outputs/ec-profiles.48.txt b/tests/07-DTT/outputs/ec-profiles.48.txt index 85d84e2..8833713 100644 --- a/tests/07-DTT/outputs/ec-profiles.48.txt +++ b/tests/07-DTT/outputs/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 1 1 12000 1.000 0.2000E-1 # COD iwarm ilarm imx ieccd: 2 5 -20 3 # COD ipec nrho istpr istpl: 1 901 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 1.11111111E-003 6.94284740E-004 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 2.22222222E-003 1.38856948E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 2 diff --git a/tests/08-STEP/outputs/ec-profiles.48.txt b/tests/08-STEP/outputs/ec-profiles.48.txt index 2cc79b5..2f2ac1f 100644 --- a/tests/08-STEP/outputs/ec-profiles.48.txt +++ b/tests/08-STEP/outputs/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 1 1 12000 1.000 0.2000 # COD iwarm ilarm imx ieccd: 2 5 -20 3 # COD ipec nrho istpr istpl: 1 501 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 2.00000000E-003 1.00806516E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 4.00000000E-003 2.01613032E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 diff --git a/tests/09-DTT-rampup/outputs/ec-profiles.48.txt b/tests/09-DTT-rampup/outputs/ec-profiles.48.txt index b85ad0f..b524db6 100644 --- a/tests/09-DTT-rampup/outputs/ec-profiles.48.txt +++ b/tests/09-DTT-rampup/outputs/ec-profiles.48.txt @@ -19,7 +19,7 @@ # COD nrayr nrayth nstep rwmax dst: 1 1 4000 1.000 0.1000 # COD iwarm ilarm imx ieccd: 2 4 20 3 # COD ipec nrho istpr istpl: 1 251 5 5 -# ρ_p ρ_t J_φ J_cdb dPdV I_cd_inside P_inside index_rt +# ρ_p ρ_t J_φ J_cd dPdV I_cd_inside P_inside index_rt 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 4.00000000E-003 3.11657387E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 8.00000000E-003 6.23314775E-003 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 0.00000000E+000 1 diff --git a/tests/12-sum-profiles/__init__.py b/tests/12-sum-profiles/__init__.py index d2f921a..dc8df52 100644 --- a/tests/12-sum-profiles/__init__.py +++ b/tests/12-sum-profiles/__init__.py @@ -82,15 +82,45 @@ class Test(TestCase): for val in ref.dtype.names: with self.subTest(value=val): - self.assertAlmostEqual(ref[val], cand[val], prec[val], - msg=f"{val} changed)") + self.assertAlmostEqual(ref[val][0], cand[val][0], + prec[val], msg=f"{val} changed)") - def test_ec_profiles(self): + def test_eccd_profiles(self): ''' - Comparing the EC radial profiles + Comparing the ECCD radial profiles ''' + from scipy.stats import wasserstein_distance as emd + import numpy as np ref = load_table(self.reference / 'sum-ec-profiles.txt') cand = load_table(self.candidate / 'sum-ec-profiles.txt') - # todo + if options.visual: + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(3, 1, sharex=True) + plt.suptitle(self.__module__ + '.test_ec_profiles') + + axes[0].set_ylabel('$J_\\text{cd}$') + axes[0].plot(ref['ρ_t'], ref['J_cd'], c='xkcd:red', + label='reference') + axes[0].plot(cand['ρ_t'], cand['J_cd'], c='xkcd:green', + ls='-.', label='candidate') + axes[0].legend() + + axes[1].set_ylabel('$dP/dV$') + axes[1].plot(ref['ρ_t'], ref['dPdV'], c='xkcd:red') + axes[1].plot(cand['ρ_t'], cand['dPdV'], c='xkcd:green', ls='-.') + + axes[2].set_xlabel('$ρ_t$') + axes[2].set_ylabel('$J_φ$') + axes[2].plot(ref['ρ_t'], ref['J_φ'], c='xkcd:red') + axes[2].plot(cand['ρ_t'], cand['J_φ'], c='xkcd:green', ls='-.') + plt.show() + + for val in ['J_cd', 'dPdV', 'J_φ']: + with self.subTest(profile=val): + y1 = abs(ref[val]) / np.sum(abs(ref[val])) + y2 = abs(cand[val]) / np.sum(abs(cand[val])) + dist = emd(ref['ρ_t'], cand['ρ_t'], y1, y2) + self.assertLess(dist, 0.01, f'{val} profile changed') diff --git a/tests/__init__.py b/tests/__init__.py index 27e7da5..46f9381 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -3,6 +3,7 @@ from typing import Any from unittest import TestCase import numpy as np +import matplotlib.pyplot as plt import shutil @@ -100,6 +101,65 @@ class GrayTest: ref[val][i], cand[val][i], prec[val], msg=f"{val} changed (ray {int(ray)})") + def test_eccd_profiles(self): + ''' + Comparing the ECCD radial profiles + ''' + from scipy.stats import wasserstein_distance as emd + import numpy as np + + try: + ref = load_table(self.reference / 'ec-profiles.48.txt') + cand = load_table(self.candidate / 'ec-profiles.48.txt') + except FileNotFoundError: + raise unittest.SkipTest("ECCD profiles not available") + + beams = np.unique(ref['index_rt']) + for index_rt, val in itertools.product(beams, ['J_cd', 'dPdV', 'J_φ']): + ref_beam = ref[ref['index_rt'] == index_rt] + cand_beam = cand[cand['index_rt'] == index_rt] + + # skip if both empty + if np.all(ref_beam[val] == 0) and np.all(cand_beam[val] == 0): + continue + + # compare with the earth mover's distance + with self.subTest(profile=val, beam=index_rt): + y1 = abs(ref_beam[val]) / np.sum(abs(ref_beam[val])) + y2 = abs(cand_beam[val]) / np.sum(abs(cand_beam[val])) + dist = emd(ref_beam['ρ_t'], cand_beam['ρ_t'], y1, y2) + self.assertLess(dist, 0.001, f'{val} profile changed') + + if options.visual: + for index_rt in beams: + ref_beam = ref[ref['index_rt'] == index_rt] + cand_beam = cand[cand['index_rt'] == index_rt] + + fig, axes = plt.subplots(3, 1, sharex=True) + fig.suptitle(self.__module__ + '.test_ec_profiles') + + axes[0].set_title(f'beam {int(index_rt)}', loc='right') + axes[0].set_ylabel('$J_\\text{cd}$') + axes[0].plot(ref_beam['ρ_t'], ref_beam['J_cd'], + c='xkcd:red', label='reference') + axes[0].plot(cand_beam['ρ_t'], cand_beam['J_cd'], + c='xkcd:green', ls='-.', label='candidate') + axes[0].legend() + + axes[1].set_ylabel('$dP/dV$') + axes[1].plot(ref_beam['ρ_t'], ref_beam['dPdV'], + c='xkcd:red') + axes[1].plot(cand_beam['ρ_t'], cand_beam['dPdV'], + c='xkcd:green', ls='-.') + + axes[2].set_xlabel('$ρ_t$') + axes[2].set_ylabel('$J_φ$') + axes[2].plot(ref_beam['ρ_t'], ref_beam['J_φ'], + c='xkcd:red') + axes[2].plot(cand_beam['ρ_t'], cand_beam['J_φ'], + c='xkcd:green', ls='-.') + plt.show() + def test_final_position(self): ''' Comparing the final position of the central ray