From 11188a7198f0c9666321afe784cf7a6b9b2b8c5e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 02:20:13 +0200 Subject: [PATCH 01/15] Update path and add CDBManager to help to export data to TAIGA --- preproc/init_from_efit.py | 53 ++++++++++++++++++++++++++------------- 1 file changed, 36 insertions(+), 17 deletions(-) diff --git a/preproc/init_from_efit.py b/preproc/init_from_efit.py index 2ec31f0..d6d18b4 100644 --- a/preproc/init_from_efit.py +++ b/preproc/init_from_efit.py @@ -39,8 +39,9 @@ def open_efit_file(self): raise def init_efit_file(self): - self.set_efit_file(get_home_directory() + '/input/cdb/' + str(self.shot_number) + - '/EFITXX/EFITXX.' + str(self.reconstruction_id) + '.h5') + self.set_efit_file(os.path.join(get_home_directory() + 'input', 'cdb', + str(self.shot_number), + 'EFITXX', 'EFITXX.' + str(self.reconstruction_id) + '.h5')) class EFITManager(EFITDataReader): @@ -74,7 +75,7 @@ def __init__(self, shot_number, time): self.efit = EFITManager(self.shot_number, self.time) self.R = self.efit.get_time_sliced_data('/output/profiles2D/r') - self.one_over_R = numpy.diag(1/self.R) + self.one_over_R = numpy.diag(1 / self.R) self.Z = self.efit.get_time_sliced_data('/output/profiles2D/z') self.poloidal_flux = [] @@ -131,25 +132,43 @@ def set_B_phi(self): class ParseToTaiga: - def __init__(self, cdb): - export_dir = get_home_directory() + '/input/fieldSpl/' + str(cdb.shot_number) + '_' + str(cdb.time) + def __init__(self, cdb, export_dir): print('Save B-spline coefficients to: ' + export_dir) + numpy.savetxt(os.path.join(export_dir, 'z.bspl'), cdb.B_R.tck[0]) + numpy.savetxt(os.path.join(export_dir, 'r.bspl'), cdb.B_R.tck[1]) + numpy.savetxt(os.path.join(export_dir, 'brad.bspl'), cdb.B_R.tck[2]) + numpy.savetxt(os.path.join(export_dir, 'bz.bspl'), cdb.B_Z.tck[2]) + numpy.savetxt(os.path.join(export_dir, 'btor.bspl'), cdb.B_phi.tck[2]) + numpy.savetxt(os.path.join(export_dir, 'psi_n.bspl'), cdb.psi_n.tck[2]) + + +class CDBManager: + export_dir = None + + def __init__(self, shot_number, time, is_forced_saved=False): + self.shot_number = shot_number + self.time = time + self.is_saved = is_forced_saved + self.init_export_dir() + if self.is_saved: + cr = CDBReader(self.shot_number, self.time) + ParseToTaiga(cr) + + def set_export_dir(self): + self.export_dir = os.path.join(get_home_directory() + 'input', 'fieldSpl', + str(self.shot_number) + '_' + str(self.time)) + + def init_export_dir(self): + self.set_export_dir() try: - os.mkdir(export_dir) - print('Create directory and write data to ' + export_dir) + self.is_saved = True + os.mkdir(self.export_dir) + print('Create B-spline coefficient directory: ' + self.export_dir) except FileExistsError: - print('Write data to ' + export_dir) + print('Existing B-spline coefficient directory: ' + self.export_dir) else: pass - numpy.savetxt(export_dir + '/z.bspl', cdb.B_R.tck[0]) - numpy.savetxt(export_dir + '/r.bspl', cdb.B_R.tck[1]) - numpy.savetxt(export_dir + '/brad.bspl', cdb.B_R.tck[2]) - numpy.savetxt(export_dir + '/bz.bspl', cdb.B_Z.tck[2]) - numpy.savetxt(export_dir + '/btor.bspl', cdb.B_phi.tck[2]) - numpy.savetxt(export_dir + '/psi_n.bspl', cdb.psi_n.tck[2]) if __name__ == "__main__": - cr = CDBReader(17178, 1097) - ParseToTaiga(cr) - + CDBManager(17178, 1097) From 7a5c72dcb2c999ed3e1329007cb547d7092cc70b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 02:40:49 +0200 Subject: [PATCH 02/15] Update mkdir --- preproc/init_from_efit.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/preproc/init_from_efit.py b/preproc/init_from_efit.py index d6d18b4..6d18648 100644 --- a/preproc/init_from_efit.py +++ b/preproc/init_from_efit.py @@ -162,12 +162,10 @@ def init_export_dir(self): self.set_export_dir() try: self.is_saved = True - os.mkdir(self.export_dir) + os.makedirs(self.export_directory) print('Create B-spline coefficient directory: ' + self.export_dir) - except FileExistsError: + except OSError as error: print('Existing B-spline coefficient directory: ' + self.export_dir) - else: - pass if __name__ == "__main__": From b0f05c03b4820b29e5cab65999a0df336f6259f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 02:49:08 +0200 Subject: [PATCH 03/15] Update mkdir --- preproc/init_from_efit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preproc/init_from_efit.py b/preproc/init_from_efit.py index 6d18648..0aa7403 100644 --- a/preproc/init_from_efit.py +++ b/preproc/init_from_efit.py @@ -164,7 +164,7 @@ def init_export_dir(self): self.is_saved = True os.makedirs(self.export_directory) print('Create B-spline coefficient directory: ' + self.export_dir) - except OSError as error: + except OSError: print('Existing B-spline coefficient directory: ' + self.export_dir) From e7a8a1de8b68adcfc64f9e7417b07b8a9e839fd7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 03:13:12 +0200 Subject: [PATCH 04/15] Update import paths --- preproc/renate_od/__init__.py | 9 +- preproc/renate_od/efit.py | 2 +- preproc/renate_od/interface.py | 113 ++++++++++++++---------- preproc/renate_od/manager.py | 2 +- preproc/renate_od/mock.py | 10 +-- preproc/renate_od/profiles.py | 8 +- preproc/renate_od/shake_flux_surface.py | 17 ++-- preproc/renate_od/thomson.py | 2 +- 8 files changed, 92 insertions(+), 71 deletions(-) diff --git a/preproc/renate_od/__init__.py b/preproc/renate_od/__init__.py index a17576e..759ecee 100644 --- a/preproc/renate_od/__init__.py +++ b/preproc/renate_od/__init__.py @@ -1,5 +1,4 @@ -import manager -from manager import RenateODManager -from beamlet import set_beamlet -from efit import EFITManager -from utils import * \ No newline at end of file +from renate_od.manager import RenateODManager +from renate_od.beamlet import set_beamlet +from renate_od.efit import EFITManager +from renate_od.utils import * diff --git a/preproc/renate_od/efit.py b/preproc/renate_od/efit.py index 00ba337..e5b8a54 100644 --- a/preproc/renate_od/efit.py +++ b/preproc/renate_od/efit.py @@ -1,6 +1,6 @@ import h5py -from utils import * +from renate_od.utils import * class EFITManager: diff --git a/preproc/renate_od/interface.py b/preproc/renate_od/interface.py index 68fe827..3fc5878 100644 --- a/preproc/renate_od/interface.py +++ b/preproc/renate_od/interface.py @@ -1,58 +1,77 @@ import matplotlib.pyplot -from manager import RenateODManager -from beamlet import set_beamlet -from efit import EFITManager -from utils import * - - -def export_beamlet_profile(export_root=os.path.join(get_home_directory(), 'input', 'ionProf'), - shot_number='17178', time='1097', species='Li', energy='80'): - z = 0 - tor = 0 - beamlet_geometry = set_beamlet(z, tor) - r = RenateODManager(beamlet_geometry, shot_number, time, species, energy) - radial_coordinate, relative_attenuation = r.get_attenuation_profile() - - export_directory = os.path.join(export_root, shot_number + '_' + time, species, energy) - try: - os.makedirs(export_directory) - print('Create directory and write data to ' + export_directory) - except OSError as error: - print('Cannot create and write to ' + export_directory) - print('Save RENATE-OD ionisation profile to: ' + export_directory) - radial_coordinate.fillna(0).to_csv(os.path.join(export_directory, 'rad.dat'), index=False, header=False) - relative_attenuation.fillna(0).to_csv(os.path.join(export_directory, 'ionyeald.dat'), index=False, header=False) - plot_attenuation_profile(shot_number, time, radial_coordinate, relative_attenuation, export_directory) - - -def get_lcfs_radial_coordinate(shot_number, time, efit_reconstruction_id=1, database_directory='input/cdb', efit_subdir='EFITXX'): +from renate_od.manager import RenateODManager +from renate_od.beamlet import set_beamlet +from renate_od.efit import EFITManager +from renate_od.utils import * + + +def get_lcfs_radial_coordinate(time, shot_number, efit_reconstruction_id=1, database_directory='input/cdb', efit_subdir='EFITXX'): data_directory = os.path.join(get_home_directory(), database_directory, str(shot_number)) efit_file = os.path.join(data_directory, efit_subdir, 'EFITXX.' + str(efit_reconstruction_id) + '.h5') efit = EFITManager(efit_file, time) return efit.get_time_sliced_data('output/separatrixGeometry/rmidplaneOut') -def plot_attenuation_profile(shot_number, time, radial_coordinate, relative_attenuation, export_directory='.'): - fig, ax = matplotlib.pyplot.subplots() - fig.set_size_inches(5, 2) - ax.plot(radial_coordinate, relative_attenuation, '-', linewidth=2) - matplotlib.pyplot.minorticks_on() - matplotlib.pyplot.grid(which='major') - matplotlib.pyplot.xlabel('$R$ [m]', labelpad=-10.5, loc='right') - matplotlib.pyplot.ylabel('neutral beam attenuation') - matplotlib.pyplot.title('COMPASS #' + shot_number + ' (' + time + ' ms)') - R_LCFS = get_lcfs_radial_coordinate(shot_number, time) - matplotlib.pyplot.axvline(R_LCFS, c='red', ls='--') - matplotlib.pyplot.text(R_LCFS+0.005, 0.45, 'LCFS', c='red', fontsize=12) - matplotlib.pyplot.savefig(os.path.join(export_directory, 'attenuation.pdf')) - matplotlib.pyplot.savefig(os.path.join(export_directory, 'attenuation.svg')) - print('Attenuation profile saved to '+ export_directory) - matplotlib.pyplot.show() +class SetProfiles: + def __init__(self, shot_number, time, species, energy, is_forced_saved=False): + self.shot_number = str(int(shot_number)) + self.time = str(int(time)) + self.species = species + self.energy = str(int(energy)) + self.export_root = os.path.join(get_home_directory(), 'input', 'ionProf') + self.is_saved = is_forced_saved + self.radial_coordinate = None + self.relative_attenuation = None + self.export_directory = None + + self.init_export_directory() + if self.is_saved: + self.get_attenuation() + self.export_beamlet_profile() + self.plot_attenuation_profile() + + def set_export_directory(self): + self.export_directory = os.path.join(self.export_root, self.shot_number + '_' + self.time, self.species, self.energy) + + def init_export_directory(self): + self.set_export_directory() + try: + os.makedirs(self.export_directory) + self.is_saved = True + print('Create attenuation profile directory: ' + self.export_directory) + except OSError: + print('Existing attenuation profile directory: ' + self.export_directory) + + def get_attenuation(self): + z = 0 + tor = 0 + beamlet_geometry = set_beamlet(z, tor) + r = RenateODManager(beamlet_geometry, self.shot_number, self.time, self.species, self.energy) + self.radial_coordinate, self.relative_attenuation = r.get_attenuation_profile() + + def export_beamlet_profile(self): + print('Save RENATE-OD ionisation profile to: ' + self.export_directory) + self.radial_coordinate.fillna(0).to_csv(os.path.join(self.export_directory, 'rad.dat'), index=False, header=False) + self.relative_attenuation.fillna(0).to_csv(os.path.join(self.export_directory, 'ionyeald.dat'), index=False, header=False) + + def plot_attenuation_profile(self): + fig, ax = matplotlib.pyplot.subplots() + fig.set_size_inches(5, 2) + ax.plot(self.radial_coordinate, self.relative_attenuation, '-', linewidth=2) + matplotlib.pyplot.minorticks_on() + matplotlib.pyplot.grid(which='major') + matplotlib.pyplot.xlabel('$R$ [m]', labelpad=-10.5, loc='right') + matplotlib.pyplot.ylabel('neutral beam attenuation') + matplotlib.pyplot.title('COMPASS #' + self.shot_number + ' (' + self.time + ' ms)') + R_LCFS = get_lcfs_radial_coordinate(time=self.time, shot_number=self.shot_number) + matplotlib.pyplot.axvline(R_LCFS, c='red', ls='--') + matplotlib.pyplot.text(R_LCFS+0.005, 0.45, 'LCFS', c='red', fontsize=12) + matplotlib.pyplot.savefig(os.path.join(self.export_directory, 'attenuation.pdf')) + matplotlib.pyplot.savefig(os.path.join(self.export_directory, 'attenuation.svg')) + print('Attenuation profile saved to '+ self.export_directory) + matplotlib.pyplot.show() if __name__ == "__main__": - a_shot_number = '17178' - a_time = '1097' - a_species = 'Li' - export_beamlet_profile(shot_number=a_shot_number, time=a_time, species=a_species) + SetProfiles(17178, 1097, 'Li', 80) diff --git a/preproc/renate_od/manager.py b/preproc/renate_od/manager.py index ec63b8c..7547591 100644 --- a/preproc/renate_od/manager.py +++ b/preproc/renate_od/manager.py @@ -2,7 +2,7 @@ import lxml.etree import sys -from profiles import * +from renate_od.profiles import * sys.path.append(os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(__file__)), os.pardir, os.pardir, 'ext', 'renate_od'))) diff --git a/preproc/renate_od/mock.py b/preproc/renate_od/mock.py index de82558..5630d8a 100644 --- a/preproc/renate_od/mock.py +++ b/preproc/renate_od/mock.py @@ -1,10 +1,10 @@ import matplotlib.pyplot -from manager import RenateODManager -from profiles import Profiles -from beamlet import BeamletGeometry, set_beamlet -from interface import get_lcfs_radial_coordinate -from utils import * +from renate_od.manager import RenateODManager +from renate_od.profiles import Profiles +from renate_od.beamlet import BeamletGeometry, set_beamlet +from renate_od.interface import get_lcfs_radial_coordinate +from renate_od.utils import * def mock_beam(shot_number='17178', time='1097', diameter=5e-3, z_length=3, tor_length=3): diff --git a/preproc/renate_od/profiles.py b/preproc/renate_od/profiles.py index 79c0339..76f6575 100644 --- a/preproc/renate_od/profiles.py +++ b/preproc/renate_od/profiles.py @@ -1,7 +1,7 @@ -from beamlet import BeamletGeometry -from thomson import ThomsonProfiles -from efit import EFITManager -from utils import * +from renate_od.beamlet import BeamletGeometry +from renate_od.thomson import ThomsonProfiles +from renate_od.efit import EFITManager +from renate_od.utils import * class Profiles: diff --git a/preproc/renate_od/shake_flux_surface.py b/preproc/renate_od/shake_flux_surface.py index 77d1477..c6c4284 100644 --- a/preproc/renate_od/shake_flux_surface.py +++ b/preproc/renate_od/shake_flux_surface.py @@ -1,12 +1,12 @@ import matplotlib.pyplot import pandas -from taiga.preproc.renate_od.interface import get_lcfs_radial_coordinate -from taiga.preproc.renate_od.manager import RenateODManager -from taiga.preproc.renate_od.beamlet import set_beamlet -from taiga.preproc.renate_od.efit import EFITManager -from taiga.preproc.renate_od.profiles import * -from taiga.preproc.renate_od.utils import * +from renate_od.interface import get_lcfs_radial_coordinate +from renate_od.manager import RenateODManager +from renate_od.beamlet import set_beamlet +from renate_od.efit import EFITManager +from renate_od.profiles import * +from renate_od.utils import * class ProfilesMock(Profiles): @@ -94,6 +94,7 @@ def shake_flux_surface_silent(species, shot_number, time, energy): (beamlet_geometry.rad[1]-beamlet_geometry.rad[2]) return numpy.max(diff) + if __name__ == "__main__": a_species = 'Li' a_shot_number = '17178' @@ -115,7 +116,9 @@ def shake_flux_surface_silent(species, shot_number, time, energy): ax.plot(energies, max_diff, '-', linewidth=2) matplotlib.pyplot.savefig(export_directory+'/maxdiff.pdf') matplotlib.pyplot.savefig(export_directory+'/maxdiff.svg') -#if __name__ == "__main__": + + +# if __name__ == "__main__": # a_species = 'Li' # a_shot_number = '17178' # a_time = '1097' diff --git a/preproc/renate_od/thomson.py b/preproc/renate_od/thomson.py index 32b8b60..7f480ef 100644 --- a/preproc/renate_od/thomson.py +++ b/preproc/renate_od/thomson.py @@ -1,6 +1,6 @@ import h5py -from utils import * +from renate_od.utils import * class ThomsonProfiles: From e4ba56c223060cd1901e4610f556e6805b5f7dd3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 03:14:08 +0200 Subject: [PATCH 05/15] Fix --- preproc/init_from_efit.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/preproc/init_from_efit.py b/preproc/init_from_efit.py index 0aa7403..93a8b15 100644 --- a/preproc/init_from_efit.py +++ b/preproc/init_from_efit.py @@ -143,16 +143,15 @@ def __init__(self, cdb, export_dir): class CDBManager: - export_dir = None - def __init__(self, shot_number, time, is_forced_saved=False): self.shot_number = shot_number self.time = time self.is_saved = is_forced_saved + self.export_dir = None self.init_export_dir() if self.is_saved: cr = CDBReader(self.shot_number, self.time) - ParseToTaiga(cr) + ParseToTaiga(cr, self.export_dir) def set_export_dir(self): self.export_dir = os.path.join(get_home_directory() + 'input', 'fieldSpl', @@ -161,8 +160,8 @@ def set_export_dir(self): def init_export_dir(self): self.set_export_dir() try: + os.makedirs(self.export_dir) self.is_saved = True - os.makedirs(self.export_directory) print('Create B-spline coefficient directory: ' + self.export_dir) except OSError: print('Existing B-spline coefficient directory: ' + self.export_dir) From b6a2a4cef1dc48c7d9e0c7c70ae5d08cad177170 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 03:14:36 +0200 Subject: [PATCH 06/15] Add preprocessor --- preproc/sd.py | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 preproc/sd.py diff --git a/preproc/sd.py b/preproc/sd.py new file mode 100644 index 0000000..2407397 --- /dev/null +++ b/preproc/sd.py @@ -0,0 +1,11 @@ +from init_from_efit import CDBManager +from renate_od.interface import SetProfiles + + +if __name__ == "__main__": + a_shot_number = 17178 + a_time = 1097 + a_species = 'Na' + an_energy = 80 + CDBManager(a_shot_number, a_time) + SetProfiles(a_shot_number, a_time, a_species, an_energy) From 1aeb89664350b32d32762f382bf23ee67ead870b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 04:08:58 +0200 Subject: [PATCH 07/15] Add core runner and parameter file --- preproc/sd.py | 54 +++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 52 insertions(+), 2 deletions(-) diff --git a/preproc/sd.py b/preproc/sd.py index 2407397..276a495 100644 --- a/preproc/sd.py +++ b/preproc/sd.py @@ -1,5 +1,56 @@ +import os +from time import localtime, strftime + from init_from_efit import CDBManager from renate_od.interface import SetProfiles +from renate_od.utils import get_home_directory + + +class SD: + def __init__(self, shot_number, time, species, energy): + self.shot_number = str(int(shot_number)) + self.time = str(int(time)) + self.species = species + self.energy = str(int(energy)) + self.runnumber = strftime("%Y%m%d%H%M%S", localtime()) + self.parameter_file = os.path.join(get_home_directory(), + "parameter_" + self.runnumber + ".sh") + + self.write_parameter_file() + self.preproc() + self.run_core() + self.delete_parameter_file() + + def preproc(self): + CDBManager(self.shot_number, self.time) + SetProfiles(self.shot_number, self.time, self.species, self.energy) + + def run_core(self): + os.chdir(get_home_directory()) + command = "./taiga_renate.exe -p=" + self.parameter_file + " -r="+self.runnumber + print(command) + os.system(command) + os.system("pwd") + + def write_parameter_file(self): + f = open(self.parameter_file, "w") + f.write("##!/usr/bin/bash\n") + f.write("shotnumber='" + self.shot_number + "'\n") + f.write("time='" + self.time + "'\n") + f.write("beammatter='" + self.species + "'\n") + f.write("energy=" + self.energy + " #keV\n") + f.write("toroidal_deflection=0 #degree; + 200 V deflection\n") + f.write("diameter=5 #mm\n") + f.write("particles=1000\n") + f.write("detector='0.6846,0.253,0.0,38,0'\n") + f.write("electric_field_module=0\n") + f.write("detector_mask='final'\n") + f.write("solver='rk'\n") + f.write("secondary_ionisation=1") + f.close() + + def delete_parameter_file(self): + os.remove(self.parameter_file) if __name__ == "__main__": @@ -7,5 +58,4 @@ a_time = 1097 a_species = 'Na' an_energy = 80 - CDBManager(a_shot_number, a_time) - SetProfiles(a_shot_number, a_time, a_species, an_energy) + SD(a_shot_number, a_time, a_species, an_energy) From e81f64c96a713a7fb8599ea44ccecd3a5ba4224c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 04:23:13 +0200 Subject: [PATCH 08/15] Update path --- preproc/renate_od/__init__.py | 8 ++++---- preproc/renate_od/efit.py | 2 +- preproc/renate_od/interface.py | 8 ++++---- preproc/renate_od/manager.py | 2 +- preproc/renate_od/mock.py | 10 +++++----- preproc/renate_od/profiles.py | 8 ++++---- preproc/renate_od/shake_flux_surface.py | 12 ++++++------ preproc/renate_od/thomson.py | 2 +- preproc/sd.py => sd.py | 6 +++--- 9 files changed, 29 insertions(+), 29 deletions(-) rename preproc/sd.py => sd.py (92%) diff --git a/preproc/renate_od/__init__.py b/preproc/renate_od/__init__.py index 759ecee..43b7d28 100644 --- a/preproc/renate_od/__init__.py +++ b/preproc/renate_od/__init__.py @@ -1,4 +1,4 @@ -from renate_od.manager import RenateODManager -from renate_od.beamlet import set_beamlet -from renate_od.efit import EFITManager -from renate_od.utils import * +from preproc.renate_od.manager import RenateODManager +from preproc.renate_od.beamlet import set_beamlet +from preproc.renate_od.efit import EFITManager +from preproc.renate_od.utils import * diff --git a/preproc/renate_od/efit.py b/preproc/renate_od/efit.py index e5b8a54..a463e6e 100644 --- a/preproc/renate_od/efit.py +++ b/preproc/renate_od/efit.py @@ -1,6 +1,6 @@ import h5py -from renate_od.utils import * +from preproc.renate_od.utils import * class EFITManager: diff --git a/preproc/renate_od/interface.py b/preproc/renate_od/interface.py index 3fc5878..f825b46 100644 --- a/preproc/renate_od/interface.py +++ b/preproc/renate_od/interface.py @@ -1,9 +1,9 @@ import matplotlib.pyplot -from renate_od.manager import RenateODManager -from renate_od.beamlet import set_beamlet -from renate_od.efit import EFITManager -from renate_od.utils import * +from preproc.renate_od.manager import RenateODManager +from preproc.renate_od.beamlet import set_beamlet +from preproc.renate_od.efit import EFITManager +from preproc.renate_od.utils import * def get_lcfs_radial_coordinate(time, shot_number, efit_reconstruction_id=1, database_directory='input/cdb', efit_subdir='EFITXX'): diff --git a/preproc/renate_od/manager.py b/preproc/renate_od/manager.py index 7547591..c68be66 100644 --- a/preproc/renate_od/manager.py +++ b/preproc/renate_od/manager.py @@ -2,7 +2,7 @@ import lxml.etree import sys -from renate_od.profiles import * +from preproc.renate_od.profiles import * sys.path.append(os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(__file__)), os.pardir, os.pardir, 'ext', 'renate_od'))) diff --git a/preproc/renate_od/mock.py b/preproc/renate_od/mock.py index 5630d8a..264f207 100644 --- a/preproc/renate_od/mock.py +++ b/preproc/renate_od/mock.py @@ -1,10 +1,10 @@ import matplotlib.pyplot -from renate_od.manager import RenateODManager -from renate_od.profiles import Profiles -from renate_od.beamlet import BeamletGeometry, set_beamlet -from renate_od.interface import get_lcfs_radial_coordinate -from renate_od.utils import * +from preproc.renate_od.manager import RenateODManager +from preproc.renate_od.profiles import Profiles +from preproc.renate_od.beamlet import BeamletGeometry, set_beamlet +from preproc.renate_od.interface import get_lcfs_radial_coordinate +from preproc.renate_od.utils import * def mock_beam(shot_number='17178', time='1097', diameter=5e-3, z_length=3, tor_length=3): diff --git a/preproc/renate_od/profiles.py b/preproc/renate_od/profiles.py index 76f6575..aacbb7b 100644 --- a/preproc/renate_od/profiles.py +++ b/preproc/renate_od/profiles.py @@ -1,7 +1,7 @@ -from renate_od.beamlet import BeamletGeometry -from renate_od.thomson import ThomsonProfiles -from renate_od.efit import EFITManager -from renate_od.utils import * +from preproc.renate_od.beamlet import BeamletGeometry +from preproc.renate_od.thomson import ThomsonProfiles +from preproc.renate_od.efit import EFITManager +from preproc.renate_od.utils import * class Profiles: diff --git a/preproc/renate_od/shake_flux_surface.py b/preproc/renate_od/shake_flux_surface.py index c6c4284..88ac506 100644 --- a/preproc/renate_od/shake_flux_surface.py +++ b/preproc/renate_od/shake_flux_surface.py @@ -1,12 +1,12 @@ import matplotlib.pyplot import pandas -from renate_od.interface import get_lcfs_radial_coordinate -from renate_od.manager import RenateODManager -from renate_od.beamlet import set_beamlet -from renate_od.efit import EFITManager -from renate_od.profiles import * -from renate_od.utils import * +from preproc.renate_od.interface import get_lcfs_radial_coordinate +from preproc.renate_od.manager import RenateODManager +from preproc.renate_od.beamlet import set_beamlet +from preproc.renate_od.efit import EFITManager +from preproc.renate_od.profiles import * +from preproc.renate_od.utils import * class ProfilesMock(Profiles): diff --git a/preproc/renate_od/thomson.py b/preproc/renate_od/thomson.py index 7f480ef..2168117 100644 --- a/preproc/renate_od/thomson.py +++ b/preproc/renate_od/thomson.py @@ -1,6 +1,6 @@ import h5py -from renate_od.utils import * +from preproc.renate_od.utils import * class ThomsonProfiles: diff --git a/preproc/sd.py b/sd.py similarity index 92% rename from preproc/sd.py rename to sd.py index 276a495..266010f 100644 --- a/preproc/sd.py +++ b/sd.py @@ -1,9 +1,9 @@ import os from time import localtime, strftime -from init_from_efit import CDBManager -from renate_od.interface import SetProfiles -from renate_od.utils import get_home_directory +from preproc.init_from_efit import CDBManager +from preproc.renate_od.interface import SetProfiles +from preproc.renate_od.utils import get_home_directory class SD: From ce429a5ed890873f8667bb11eaeab4d147e2437d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 04:40:50 +0200 Subject: [PATCH 09/15] Add plotter --- plotter/detector.py | 113 +++++++++----------- plotter/detector_plane.py | 211 +++++++++++++++++--------------------- sd.py | 16 ++- 3 files changed, 155 insertions(+), 185 deletions(-) diff --git a/plotter/detector.py b/plotter/detector.py index d5e7ab6..a3b08fb 100644 --- a/plotter/detector.py +++ b/plotter/detector.py @@ -5,72 +5,57 @@ import matplotlib.cm as colormap import scipy.stats as st -try: - shotnumber = sys.argv[1]; -except: - print 'Shot number is not given as proper input parameter. Default value (17178) is used.' - shotnumber = '17178' -try: - time = sys.argv[2]; -except: - print 'Time is not given as proper input parameter. Default value (1005) is used.' - time = '1005' - -try: - runnumber = sys.argv[3]; -except: - print 'Runnumber is not given as proper input parameter. Default value (402) is used.' - runnumber = '402' - -results_folder = 'results/'+shotnumber+'_'+time+'/'+runnumber+'/' - -try: - C = np.loadtxt(results_folder+'detector/cellcounter.dat') - X = np.loadtxt(results_folder+'detector/detx') - Y = np.loadtxt(results_folder+'detector/dety') -except: - print 'Invalid input folder: '+results_folder - - -xmin = X[0,0] -xmax = X[-1,-1] -ymin = Y[0,0] -ymax = Y[-1,-1] -xmean = (X[:,0]+X[:,1])/2 -ymean = (Y[:,0]+Y[:,1])/2 -cmax = C.max() - -try: - fig = pl.figure() - ax = fig.gca() - ax.set_xlim(xmin, xmax) - ax.set_ylim(ymin, ymax) - ax.set_aspect('equal') - - # Contourf plot - pl.contourf(xmean, ymean, C, cmap='afmhot', corner_mask=False) - - # Detector cell matrix plot - for i in range (X[:,0].size): - for j in range (Y[:,0].size): - rect = patches.Rectangle((X[i,0], Y[j,0]), X[i,1]-X[i,0], Y[j,1]-Y[j,0], linewidth=1, - edgecolor='k', facecolor=colormap.afmhot(C[j,i]/cmax)) - ax.add_patch(rect) - - print(C) - - pl.xlabel(r"$x \mathrm{ [mm]}$") - pl.ylabel(r"$y \mathrm{ [mm]}$") - pl.title(r"COMPASS #"+shotnumber+" (t="+time+" s)") +def detector(shotnumber, time, runnumber): + results_folder = 'results/'+shotnumber+'_'+time+'/'+runnumber+'/' + + try: + C = np.loadtxt(results_folder+'detector/cellcounter.dat') + X = np.loadtxt(results_folder+'detector/detx') + Y = np.loadtxt(results_folder+'detector/dety') + except: + print ('Invalid input folder: '+results_folder) + + + xmin = X[0,0] + xmax = X[-1,-1] + ymin = Y[0,0] + ymax = Y[-1,-1] + xmean = (X[:,0]+X[:,1])/2 + ymean = (Y[:,0]+Y[:,1])/2 + cmax = C.max() + try: - print 'Save plot to '+results_folder+'detpy2_'+shotnumber+'_'+time+'.pdf' - pl.savefig(results_folder+'detpy2_'+shotnumber+'_'+time+'.pdf') - pl.clf() + fig = pl.figure() + ax = fig.gca() + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + ax.set_aspect('equal') + + # Contourf plot + pl.contourf(xmean, ymean, C, cmap='afmhot', corner_mask=False) + + # Detector cell matrix plot + for i in range (X[:,0].size): + for j in range (Y[:,0].size): + rect = patches.Rectangle((X[i,0], Y[j,0]), X[i,1]-X[i,0], Y[j,1]-Y[j,0], linewidth=1, + edgecolor='k', facecolor=colormap.afmhot(C[j,i]/cmax)) + ax.add_patch(rect) + + print(C) + + pl.xlabel(r"$x \mathrm{ [mm]}$") + pl.ylabel(r"$y \mathrm{ [mm]}$") + pl.title(r"COMPASS #"+shotnumber+" (t="+time+" s)") + + try: + print('Save plot to '+results_folder+'detpy2_'+shotnumber+'_'+time+'.pdf') + pl.savefig(results_folder+'detpy2_'+shotnumber+'_'+time+'.pdf') + pl.clf() + except: + print('Unable to save to : '+results_folder) + pl.show() except: - print 'Unable to save to : '+results_folder - pl.show() -except: - print 'Unable to plot: '+results_folder + print('Unable to plot: '+results_folder) diff --git a/plotter/detector_plane.py b/plotter/detector_plane.py index 49e1574..52cf1de 100644 --- a/plotter/detector_plane.py +++ b/plotter/detector_plane.py @@ -3,120 +3,99 @@ import matplotlib.pyplot as pl import scipy.stats as st -max_distance = 1e-4 - -try: - shotnumber = sys.argv[1]; -except: - print 'Shot number is not given as proper input parameter. Default value (17178) is used.' - shotnumber = '17178' - -try: - time = sys.argv[2]; -except: - print 'Time is not given as proper input parameter. Default value (1005) is used.' - time = '1005' - -try: - runnumber = sys.argv[3]; -except: - print 'Runnumber is not given as proper input parameter. Default value (402) is used.' - runnumber = '402' - -try: - detector_par = sys.argv[4]; -except: - print 'Detector parameters is not given as proper input parameter. Default value (0.685,0.22,0,38,0) is used.' - detector_par = '0.685,0.22,0,38,0' - -detector_par_list = detector_par.split(',') - -try: - R0 = float(detector_par_list[0]); -except: - R0 = 0.685 - -try: - Z0 = float(detector_par_list[1]); -except: - Z0 = 0.23 - -try: - T0 = float(detector_par_list[2]); -except: - T0 = 0.0 - -try: - detector_RZ_angle = float(detector_par_list[3]); -except: - detector_RZ_angle = 38.0 - -try: - detector_ZT_angle = float(detector_par_list[4]); -except: - detector_ZT_angle = 0.0 - -results_folder = 'results/'+shotnumber+'_'+time+'/'+runnumber+'/' - -try: - R = np.loadtxt(results_folder+'rad.dat') #, delimiter="\t" - Z = np.loadtxt(results_folder+'z.dat') - T = np.loadtxt(results_folder+'tor.dat') -except: - print 'Invalid input folder: '+results_folder - -try: - particle_on_detector = np.abs((Z-Z0) *np.tan((detector_RZ_angle/180.0)*np.pi) + (R-R0)) < max_distance - - x = T[particle_on_detector] - y = Z[particle_on_detector] + +def detector_plane(shotnumber, time, runnumber, detector_par): + + max_distance = 1e-4 + + detector_par_list = detector_par.split(',') + + try: + R0 = float(detector_par_list[0]); + except: + R0 = 0.685 + + try: + Z0 = float(detector_par_list[1]); + except: + Z0 = 0.23 + + try: + T0 = float(detector_par_list[2]); + except: + T0 = 0.0 + + try: + detector_RZ_angle = float(detector_par_list[3]); + except: + detector_RZ_angle = 38.0 + + try: + detector_ZT_angle = float(detector_par_list[4]); + except: + detector_ZT_angle = 0.0 + + results_folder = 'results/'+shotnumber+'_'+time+'/'+runnumber+'/' + + try: + R = np.loadtxt(results_folder+'rad.dat') #, delimiter="\t" + Z = np.loadtxt(results_folder+'z.dat') + T = np.loadtxt(results_folder+'tor.dat') + except: + print ('Invalid input folder: '+results_folder) + + try: + particle_on_detector = np.abs((Z-Z0) *np.tan((detector_RZ_angle/180.0)*np.pi) + (R-R0)) < max_distance + + x = T[particle_on_detector] + y = Z[particle_on_detector] + xmin, xmax = -4e-2, 4e-2 + ymin, ymax = 18e-2, 26e-2 + + # Peform the kernel density estimate + xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] + positions = np.vstack([xx.ravel(), yy.ravel()]) + values = np.vstack([x, y]) + kernel = st.gaussian_kde(values,0.05) + f = np.reshape(kernel(positions).T, xx.shape) + + fig = pl.figure() + ax = fig.gca() + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + ax.set_aspect('equal') + + # Contourf plot + #cfset = ax.contourf(xx, yy, f, 500, cmap='hot') # Blues, afmhot; magma, inferno + ax.hist2d(x, y, bins=250, range = [[xmin,xmax],[ymin,ymax]], cmap='afmhot') + # Contour plot + # Label plot + + pl.xlabel(r"$T \mathrm{ [m]}$") + pl.ylabel(r"$Z \mathrm{ [m]}$") + except: + print ('Unable to plot: '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') + + try: + print ('Save plot to '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') + pl.savefig(results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') + pl.clf() + + except: + print ('Unable to save to : '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') + pl.show() + xmin, xmax = -4e-2, 4e-2 - ymin, ymax = 18e-2, 26e-2 - - # Peform the kernel density estimate - xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] - positions = np.vstack([xx.ravel(), yy.ravel()]) - values = np.vstack([x, y]) - kernel = st.gaussian_kde(values,0.05) - f = np.reshape(kernel(positions).T, xx.shape) - - fig = pl.figure() - ax = fig.gca() - ax.set_xlim(xmin, xmax) - ax.set_ylim(ymin, ymax) - ax.set_aspect('equal') - - # Contourf plot - #cfset = ax.contourf(xx, yy, f, 500, cmap='hot') # Blues, afmhot; magma, inferno - ax.hist2d(x, y, bins=250, range = [[xmin,xmax],[ymin,ymax]], cmap='afmhot') - # Contour plot - # Label plot - - pl.xlabel(r"$T \mathrm{ [m]}$") - pl.ylabel(r"$Z \mathrm{ [m]}$") -except: - print 'Unable to plot: '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf' - -try: - print 'Save plot to '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf' - pl.savefig(results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') - pl.clf() - -except: - print 'Unable to save to : '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf' - pl.show() - -xmin, xmax = -4e-2, 4e-2 -ymin, ymax = -4e-2, 4e-2 - -#x = T[particle_on_detector]-T0 -#y = (Z[particle_on_detector]-Z0) / np.sin( (detector_RZ_angle/180.0)*np.pi ) -#nbins=250 -#fig = pl.figure() -#ax = fig.gca() -#ax.set_xlim(xmin, xmax) -#ax.set_ylim(ymin, ymax) -#ax.set_aspect('equal') -#cfset = ax.hist2d(x, y, bins=nbins, range = [[xmin,xmax],[ymin,ymax]], cmap='afmhot') -#pl.colorbar(cfset[3]),ax=ax) -#pl.show() + ymin, ymax = -4e-2, 4e-2 + + #x = T[particle_on_detector]-T0 + #y = (Z[particle_on_detector]-Z0) / np.sin( (detector_RZ_angle/180.0)*np.pi ) + #nbins=250 + #fig = pl.figure() + #ax = fig.gca() + #ax.set_xlim(xmin, xmax) + #ax.set_ylim(ymin, ymax) + #ax.set_aspect('equal') + #cfset = ax.hist2d(x, y, bins=nbins, range = [[xmin,xmax],[ymin,ymax]], cmap='afmhot') + #pl.colorbar(cfset[3]),ax=ax) + #pl.show() diff --git a/sd.py b/sd.py index 266010f..4ba3ed9 100644 --- a/sd.py +++ b/sd.py @@ -4,14 +4,17 @@ from preproc.init_from_efit import CDBManager from preproc.renate_od.interface import SetProfiles from preproc.renate_od.utils import get_home_directory +from plotter.detector import detector +from plotter.detector_plane import detector_plane class SD: - def __init__(self, shot_number, time, species, energy): + def __init__(self, shot_number, time, species, energy, detector_par): self.shot_number = str(int(shot_number)) self.time = str(int(time)) self.species = species self.energy = str(int(energy)) + self.detector_par = detector_par self.runnumber = strftime("%Y%m%d%H%M%S", localtime()) self.parameter_file = os.path.join(get_home_directory(), "parameter_" + self.runnumber + ".sh") @@ -20,6 +23,8 @@ def __init__(self, shot_number, time, species, energy): self.preproc() self.run_core() self.delete_parameter_file() + detector(self.shot_number, self.time, self.runnumber) + detector_plane(self.shot_number, self.time, self.runnumber, self.detector_par) def preproc(self): CDBManager(self.shot_number, self.time) @@ -42,7 +47,7 @@ def write_parameter_file(self): f.write("toroidal_deflection=0 #degree; + 200 V deflection\n") f.write("diameter=5 #mm\n") f.write("particles=1000\n") - f.write("detector='0.6846,0.253,0.0,38,0'\n") + f.write("detector='" + self.detector_par + "'\n") f.write("electric_field_module=0\n") f.write("detector_mask='final'\n") f.write("solver='rk'\n") @@ -56,6 +61,7 @@ def delete_parameter_file(self): if __name__ == "__main__": a_shot_number = 17178 a_time = 1097 - a_species = 'Na' - an_energy = 80 - SD(a_shot_number, a_time, a_species, an_energy) + a_species = 'Li' + an_energy = 70 + a_detector = "0.6846,0.253,0.0,38,0" + SD(a_shot_number, a_time, a_species, an_energy, a_detector) From e5b1a67e600351cd6f50d47977c8e07b733bfa78 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 23:22:47 +0200 Subject: [PATCH 10/15] Fix plotters --- plotter/detector.py | 10 +- plotter/detector_plane.py | 8 +- plotter/traj_plotter.py | 285 +++++++++++++++++++------------------- sd.py | 47 +++++-- 4 files changed, 187 insertions(+), 163 deletions(-) diff --git a/plotter/detector.py b/plotter/detector.py index a3b08fb..d53c4cf 100644 --- a/plotter/detector.py +++ b/plotter/detector.py @@ -6,7 +6,7 @@ import scipy.stats as st -def detector(shotnumber, time, runnumber): +def detector(shotnumber, time, runnumber, is_debug=False): results_folder = 'results/'+shotnumber+'_'+time+'/'+runnumber+'/' @@ -25,6 +25,9 @@ def detector(shotnumber, time, runnumber): xmean = (X[:,0]+X[:,1])/2 ymean = (Y[:,0]+Y[:,1])/2 cmax = C.max() + + if cmax == 0: + cmax = 1 try: fig = pl.figure() @@ -42,8 +45,9 @@ def detector(shotnumber, time, runnumber): rect = patches.Rectangle((X[i,0], Y[j,0]), X[i,1]-X[i,0], Y[j,1]-Y[j,0], linewidth=1, edgecolor='k', facecolor=colormap.afmhot(C[j,i]/cmax)) ax.add_patch(rect) - - print(C) + + if is_debug: + print(C) pl.xlabel(r"$x \mathrm{ [mm]}$") pl.ylabel(r"$y \mathrm{ [mm]}$") diff --git a/plotter/detector_plane.py b/plotter/detector_plane.py index 52cf1de..34ddf7b 100644 --- a/plotter/detector_plane.py +++ b/plotter/detector_plane.py @@ -42,7 +42,7 @@ def detector_plane(shotnumber, time, runnumber, detector_par): Z = np.loadtxt(results_folder+'z.dat') T = np.loadtxt(results_folder+'tor.dat') except: - print ('Invalid input folder: '+results_folder) + print('Invalid input folder: ' + results_folder) try: particle_on_detector = np.abs((Z-Z0) *np.tan((detector_RZ_angle/180.0)*np.pi) + (R-R0)) < max_distance @@ -74,15 +74,15 @@ def detector_plane(shotnumber, time, runnumber, detector_par): pl.xlabel(r"$T \mathrm{ [m]}$") pl.ylabel(r"$Z \mathrm{ [m]}$") except: - print ('Unable to plot: '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') + print('Unable to plot: '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') try: - print ('Save plot to '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') + print('Save plot to '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') pl.savefig(results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') pl.clf() except: - print ('Unable to save to : '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') + print('Unable to save to : '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') pl.show() xmin, xmax = -4e-2, 4e-2 diff --git a/plotter/traj_plotter.py b/plotter/traj_plotter.py index 5984a08..c5a2a24 100644 --- a/plotter/traj_plotter.py +++ b/plotter/traj_plotter.py @@ -5,166 +5,165 @@ import h5py import matplotlib.pyplot as plt -numberOfCmdArgs = len(sys.argv) - -rootFolder = os.path.abspath(os.getcwd()) -print(rootFolder) -resultFolder = rootFolder+'/results'; -fieldFolder = rootFolder+'/input/fieldGrid'; -cdbFolder = rootFolder+'/input/cdb'; -ionProfileFolder =rootFolder+'/input/ionProf'; - -print(numberOfCmdArgs) -if (numberOfCmdArgs > 1): - shotAndTime = sys.argv[1] -else: - shotAndTime = '11774_1000' - -if (numberOfCmdArgs > 2): - runnumber = sys.argv[2] -else: - runnumber = '0'; - -workingFolder = resultFolder+'/'+shotAndTime+'/'+runnumber - -#ionisation -ionR = numpy.genfromtxt(fname=ionProfileFolder+'/'+shotAndTime+'/rad.dat') -ionY = numpy.genfromtxt(fname=ionProfileFolder+'/'+shotAndTime+'/ionyeald.dat') -ionX1 = numpy.gradient(ionY); -ionF = interp1d(ionR, -ionX1/numpy.amax(numpy.absolute(ionX1))); - -plt.plot(ionR, ionF(ionR)) - -try: - print ('Save plot to '+workingFolder+'/ion_'+shotAndTime+'.pdf') - plt.savefig(workingFolder+'/ion_'+shotAndTime+'.pdf') - plt.clf() -except: - print ('Unable to save to : '+workingFolder+'/ion_'+shotAndTime+'.pdf') - plt.show() +def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy): + rootFolder = os.path.abspath(os.getcwd()) + resultFolder = rootFolder+'/results' + fieldFolder = rootFolder+'/input/fieldGrid' + cdbFolder = rootFolder+'/input/cdb' + ionProfileFolder =rootFolder+'/input/ionProf' + + shotAndTime = shotnumber + '_' + time + + workingFolder = resultFolder+'/'+shotAndTime+'/'+runnumber + + #ionisation + ionR = numpy.genfromtxt(fname=ionProfileFolder+'/'+shotAndTime+'/'+species+'/'+energy+'/rad.dat') + ionY = numpy.genfromtxt(fname=ionProfileFolder+'/'+shotAndTime+'/'+species+'/'+energy+'/ionyeald.dat') + ionX1 = numpy.gradient(ionY) + ionF = interp1d(ionR, -ionX1/numpy.amax(numpy.absolute(ionX1))) + + plt.plot(ionR, ionF(ionR)) -# trajectories -print(resultFolder+'/'+shotAndTime+'/'+runnumber+'/t_rad.dat') -t_rad = numpy.genfromtxt(fname=workingFolder+'/t_rad.dat') -t_z = numpy.genfromtxt(fname=workingFolder+'/t_z.dat') -t_tor = numpy.genfromtxt(fname=workingFolder+'/t_tor.dat') + try: + print('Save plot to '+workingFolder+'/ion_'+shotAndTime+'.pdf') + plt.savefig(workingFolder+'/ion_'+shotAndTime+'.pdf') + plt.clf() -numberOfIons = numpy.size(t_rad, 1) + except: + print('Unable to save to : '+workingFolder+'/ion_'+shotAndTime+'.pdf') + plt.show() -print(fieldFolder+'/'+shotAndTime+'/rcord.dat') -rcord = numpy.genfromtxt(fname=fieldFolder+'/'+shotAndTime+'/rcord.dat') -zcord = numpy.genfromtxt(fname=fieldFolder+'/'+shotAndTime+'/zcord.dat') -psi2 = numpy.genfromtxt(fname=fieldFolder+'/'+shotAndTime+'/psi2.dat') + # trajectories + t_rad = numpy.genfromtxt(fname=workingFolder+'/t_rad.dat') + t_z = numpy.genfromtxt(fname=workingFolder+'/t_z.dat') + t_tor = numpy.genfromtxt(fname=workingFolder+'/t_tor.dat') + cellid = numpy.genfromtxt(fname=workingFolder+'/detector/cellid.dat') -shotnumber = shotAndTime.split('_')[0] -timeSlice = shotAndTime.split('_')[1] + rcord = numpy.genfromtxt(fname=fieldFolder+'/'+shotAndTime+'/rcord.dat') + zcord = numpy.genfromtxt(fname=fieldFolder+'/'+shotAndTime+'/zcord.dat') + psi2 = numpy.genfromtxt(fname=fieldFolder+'/'+shotAndTime+'/psi2.dat') -hdf5Path = cdbFolder+'/'+shotnumber+'/EFITXX/EFITXX.1.h5' -print(hdf5Path) -hdf5File = h5py.File(hdf5Path, 'r') -hdf5Times = hdf5File.get('time').value -hdf5PsiBoundary = hdf5File.get('output/globalParameters/psiBoundary').value -hdf5PsiAxis = hdf5File.get('output/globalParameters/psiAxis').value -hdf5BoundaryR = hdf5File.get('output/separatrixGeometry/boundaryCoordsR').value -hdf5BoundaryZ = hdf5File.get('output/separatrixGeometry/boundaryCoordsZ').value + hdf5Path = cdbFolder+'/'+shotnumber+'/EFITXX/EFITXX.1.h5' + hdf5File = h5py.File(hdf5Path, 'r') + hdf5Times = hdf5File.get('time')[()] + hdf5PsiBoundary = hdf5File.get('output/globalParameters/psiBoundary')[()] + hdf5PsiAxis = hdf5File.get('output/globalParameters/psiAxis')[()] + hdf5BoundaryR = hdf5File.get('output/separatrixGeometry/boundaryCoordsR')[()] + hdf5BoundaryZ = hdf5File.get('output/separatrixGeometry/boundaryCoordsZ')[()] -hdf5Index = int(numpy.argwhere(hdf5Times == float(timeSlice)/1000.)) + hdf5Index = int(numpy.argwhere(hdf5Times == float(time)/1000.)) -psiAxis = hdf5PsiAxis[hdf5Index] -psiBoundary = hdf5PsiBoundary[hdf5Index] -boundaryR = hdf5BoundaryR[hdf5Index,:] -boundaryZ = hdf5BoundaryZ[hdf5Index,:] + psiAxis = hdf5PsiAxis[hdf5Index] + psiBoundary = hdf5PsiBoundary[hdf5Index] + boundaryR = hdf5BoundaryR[hdf5Index,:] + boundaryZ = hdf5BoundaryZ[hdf5Index,:] -PSI = numpy.transpose(psi2)/psiAxis + PSI = numpy.transpose(psi2)/psiAxis -R, Z = numpy.meshgrid(rcord, zcord) + R, Z = numpy.meshgrid(rcord, zcord) -if (numberOfCmdArgs > 3): - t_index = sys.argv[2] + t_index = cellid >= 0 t_rad = t_rad[:,t_index] t_z = t_z [:,t_index] t_tor = t_tor[:,t_index] + numberOfIons = numpy.size(t_rad, 1) -#ion histogram (real) -fig = plt.figure() -plt.hist(t_rad[0,:], density=True, bins=50) -plt.xlabel(r'$R$ [m]') -plt.title (r'COMPASS #'+shotnumber+'\n ($t = $'+timeSlice+' ms)') + #ion histogram (real) + fig = plt.figure() + plt.hist(t_rad[0,:], density=True, bins=50) + plt.xlabel(r'$R$ [m]') + plt.title (r'COMPASS #'+shotnumber+'\n ($t = $'+time+' ms)') -try: - print ('Save plot to '+workingFolder+'/start_'+shotAndTime+'.pdf') - plt.savefig(workingFolder+'/start_'+shotAndTime+'.pdf') - plt.clf() + try: + print('Save plot to '+workingFolder+'/start_'+shotAndTime+'.pdf') + plt.savefig(workingFolder+'/start_'+shotAndTime+'.pdf') + plt.clf() -except: - print ('Unable to save to : '+workingFolder+'/start_'+shotAndTime+'.pdf') - plt.show() - -#detector -plt.clf() -fig, axs = plt.subplots(2, 2) -axs[0,0].plot(t_rad[-1,:], t_z[-1,:], '.') -#axs[0,0].set_xlabel(r'$R$ [m]') -axs[0,0].set_ylabel(r'$Z$ [m]') -axs[0,0].yaxis.set_ticks_position('both') - -axs[0,1].plot(t_tor[-1,:], t_z[-1,:], '.') -axs[0,1].set_xlabel(r'$T$ [m]') -axs[0,1].set_ylabel(r'$Z$ [m]') -axs[0,1].yaxis.tick_right() -axs[0,1].yaxis.set_ticks_position('both') - -axs[1,0].plot(t_rad[-1,:], t_tor[-1,:], '.') -axs[1,0].set_xlabel(r'$R$ [m]') -axs[1,0].set_ylabel(r'$T$ [m]') -axs[1,0].yaxis.set_ticks_position('both') - -axs[1,1].axis('off') - -fig.suptitle(r'COMPASS #'+shotnumber+' ($t = $'+timeSlice+' ms)') -#plt.colorbar() - -try: - print ('Save plot to '+workingFolder+'/end_'+shotAndTime+'.pdf') - plt.savefig(workingFolder+'/end_'+shotAndTime+'.pdf') - plt.clf() + except: + print('Unable to save to : '+workingFolder+'/start_'+shotAndTime+'.pdf') + plt.show() -except: - print ('Unable to save to : '+workingFolder+'/end_'+shotAndTime+'.pdf') - plt.show() - -#2D figure -plt.clf() -fig = plt.figure() -#plt.contourf(R, Z, -PSI, levels=numpy.arange(-1.05,2,.05), cmap='Spectral') -#plt.contourf(R, Z, PSI, levels=numpy.arange(-2,2,.05), cmap='coolwarm') -#plt.contourf(R, Z, PSI, levels=numpy.arange(-2,2,.05), cmap='twilight_shifted') -plt.contourf(R, Z, PSI, levels=numpy.arange(-2,2,.1), cmap='GnBu') -plt.plot(boundaryR, boundaryZ, color=(0,0,0)) -for i in range(numberOfIons): - i1 = 1 - i2 = ionF(t_rad[1,i]) - plt.plot([rcord[-1], t_rad[1,i]], [t_z[0,i], t_z[0,i]], color=(i1, i1, 0)) - if i2<0: - i2=0 - if i2>1: - i2=1 - plt.plot(t_rad[:,i], t_z[:,i], color=(i2, 0, 0)) - -plt.xlabel(r'$R$ [m]') -plt.ylabel(r'$Z$ [m]') -plt.title (r'COMPASS #'+shotnumber+'\n ($t = $'+timeSlice+' ms)') -ax = fig.add_subplot(111) -ax.set_aspect('equal', 'box') -#plt.colorbar() - -try: - print ('Save plot to '+workingFolder+'/traj_'+shotAndTime+'.pdf') - plt.savefig(workingFolder+'/traj_'+shotAndTime+'.pdf') + #detector plt.clf() - -except: - print ('Unable to save to : '+workingFolder+'/traj_'+shotAndTime+'.pdf') - plt.show() + fig, axs = plt.subplots(2, 2) + axs[0,0].plot(t_rad[-1,:], t_z[-1,:], '.') + #axs[0,0].set_xlabel(r'$R$ [m]') + axs[0,0].set_ylabel(r'$Z$ [m]') + axs[0,0].yaxis.set_ticks_position('both') + + axs[0,1].plot(t_tor[-1,:], t_z[-1,:], '.') + axs[0,1].set_xlabel(r'$T$ [m]') + axs[0,1].set_ylabel(r'$Z$ [m]') + axs[0,1].yaxis.tick_right() + axs[0,1].yaxis.set_ticks_position('both') + + axs[1,0].plot(t_rad[-1,:], t_tor[-1,:], '.') + axs[1,0].set_xlabel(r'$R$ [m]') + axs[1,0].set_ylabel(r'$T$ [m]') + axs[1,0].yaxis.set_ticks_position('both') + + axs[1,1].axis('off') + + fig.suptitle(r'COMPASS #'+shotnumber+' ($t = $'+time+' ms)') + #plt.colorbar() + + try: + print('Save plot to '+workingFolder+'/end_'+shotAndTime+'.pdf') + plt.savefig(workingFolder+'/end_'+shotAndTime+'.pdf') + plt.clf() + + except: + print('Unable to save to : '+workingFolder+'/end_'+shotAndTime+'.pdf') + plt.show() + + #2D figure + plt.clf() + fig = plt.figure() + ax = fig.add_subplot(111) + #plt.contourf(R, Z, -PSI, levels=numpy.arange(-1.05,2,.05), cmap='Spectral') + #plt.contourf(R, Z, PSI, levels=numpy.arange(-2,2,.05), cmap='coolwarm') + #plt.contourf(R, Z, PSI, levels=numpy.arange(-2,2,.05), cmap='twilight_shifted') + plt.contourf(R, Z, PSI, levels=numpy.arange(-2,2,.1), cmap='GnBu') + plt.plot(boundaryR, boundaryZ, color=(0,0,0)) + for i in range(numberOfIons): + i1 = 1 + i2 = numpy.asscalar(ionF(t_rad[1,i])) + plt.plot([rcord[-1], t_rad[1,i]], [t_z[0,i], t_z[0,i]], color=(i1, i1, 0)) + if i2<0: + i2=0 + if i2>1: + i2=1 + plt.plot(t_rad[:,i], t_z[:,i], color=(i2, 0, 0)) + + # plot detector + try: + detector_y_size = 12.2e-3 + detector_x0, detector_y0, detector_z0, detector_phi, detector_psi = map(float, detector_par.split(',')) + + detector_x_tilt = detector_y_size * numpy.sin(detector_phi / 180 * numpy.pi) + detector_y_tilt = detector_y_size * numpy.cos(detector_phi / 180 * numpy.pi) + plt.plot([detector_x0, detector_x0], [detector_y0, max(zcord)], color=(0, 0, 1)) + plt.plot([detector_x0, detector_x0], [detector_y0, max(zcord)], color=(0, 0, 1)) + plt.plot([detector_x0 - detector_x_tilt, detector_x0 + detector_x_tilt], + [detector_y0 + detector_y_tilt, detector_y0 - detector_y_tilt], + color=(0, 0, 1), linewidth=2) + except: + print('Invalid detector coordinates.') + pass + + plt.xlabel(r'$R$ [m]') + plt.ylabel(r'$Z$ [m]') + plt.title (r'COMPASS #'+shotnumber+'\n ($t = $'+time+' ms)') + ax.set_aspect('equal', 'box') + #plt.colorbar() + + try: + print('Save plot to '+workingFolder+'/traj_'+shotAndTime+'.pdf') + plt.savefig(workingFolder+'/traj_'+shotAndTime+'.pdf') + plt.clf() + + except: + print('Unable to save to : '+workingFolder+'/traj_'+shotAndTime+'.pdf') + plt.show() diff --git a/sd.py b/sd.py index 4ba3ed9..53f1e05 100644 --- a/sd.py +++ b/sd.py @@ -6,25 +6,30 @@ from preproc.renate_od.utils import get_home_directory from plotter.detector import detector from plotter.detector_plane import detector_plane +from plotter.traj_plotter import traj_plotter class SD: - def __init__(self, shot_number, time, species, energy, detector_par): + def __init__(self, shot_number, time, species, energy, detector_par, + runnumber=None, is_trajectory_detailed=False): self.shot_number = str(int(shot_number)) self.time = str(int(time)) self.species = species self.energy = str(int(energy)) self.detector_par = detector_par - self.runnumber = strftime("%Y%m%d%H%M%S", localtime()) - self.parameter_file = os.path.join(get_home_directory(), - "parameter_" + self.runnumber + ".sh") + self.runnumber = runnumber + self.is_trajectory_detailed = is_trajectory_detailed + self.is_run_simulation = True + self.set_runnumber() + if self.is_run_simulation: + self.parameter_file = os.path.join(get_home_directory(), + "parameter_" + self.runnumber + ".sh") - self.write_parameter_file() - self.preproc() - self.run_core() - self.delete_parameter_file() - detector(self.shot_number, self.time, self.runnumber) - detector_plane(self.shot_number, self.time, self.runnumber, self.detector_par) + self.write_parameter_file() + self.preproc() + self.run_core() + self.delete_parameter_file() + self.run_plotters() def preproc(self): CDBManager(self.shot_number, self.time) @@ -32,10 +37,24 @@ def preproc(self): def run_core(self): os.chdir(get_home_directory()) - command = "./taiga_renate.exe -p=" + self.parameter_file + " -r="+self.runnumber + command = "./taiga_renate.exe -p=" + self.parameter_file + " -r=" + self.runnumber + if self.is_trajectory_detailed: + command += " --fulltrace" print(command) os.system(command) - os.system("pwd") + + def run_plotters(self): + os.chdir(get_home_directory()) + detector(self.shot_number, self.time, self.runnumber) + detector_plane(self.shot_number, self.time, self.runnumber, self.detector_par) + if self.is_trajectory_detailed: + traj_plotter(self.shot_number, self.time, self.runnumber, self.detector_par, self.species, self.energy) + + def set_runnumber(self): + if self.runnumber is None: + self.runnumber = strftime("%Y%m%d%H%M%S", localtime()) + else: + self.is_run_simulation = False def write_parameter_file(self): f = open(self.parameter_file, "w") @@ -64,4 +83,6 @@ def delete_parameter_file(self): a_species = 'Li' an_energy = 70 a_detector = "0.6846,0.253,0.0,38,0" - SD(a_shot_number, a_time, a_species, an_energy, a_detector) + a_runnumber = '20240523214835' + is_plot_trajectory = True + SD(a_shot_number, a_time, a_species, an_energy, a_detector, a_runnumber, is_trajectory_detailed=is_plot_trajectory) From 8eeb467aeb5c032104336054f559770246800a75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Thu, 23 May 2024 23:23:33 +0200 Subject: [PATCH 11/15] Remove old plotters --- plotter/calculate_larmor_radius.m | 111 ------------------------------ plotter/plot3dtaiga.m | 28 -------- plotter/traj_plotter.m | 58 ---------------- 3 files changed, 197 deletions(-) delete mode 100644 plotter/calculate_larmor_radius.m delete mode 100644 plotter/plot3dtaiga.m delete mode 100644 plotter/traj_plotter.m diff --git a/plotter/calculate_larmor_radius.m b/plotter/calculate_larmor_radius.m deleted file mode 100644 index 9a7b829..0000000 --- a/plotter/calculate_larmor_radius.m +++ /dev/null @@ -1,111 +0,0 @@ -clear all - mainfolder = '../results'; - shotnumber = 'test_1'; runnumber = '120'; - shotnumber = 'test_2'; runnumber = '128'; - shotnumber = 'test_5'; runnumber = '129'; - shotnumber = 'test_6'; runnumber = '130'; - shotnumber = 'test_10'; runnumber = '132'; - - - load([mainfolder,'/',shotnumber,'/',runnumber,'/t_rad.dat']); - load([mainfolder,'/',shotnumber,'/',runnumber,'/t_z.dat']); - load([mainfolder,'/',shotnumber,'/',runnumber,'/t_tor.dat']); - - - t_rad = t_rad'; - t_z = t_z'; - t_tor = t_tor'; - -if false - t_rad = t_rad([1,1],:); - t_z = t_z([1,1],:); - t_tor = t_tor([1,1],:); -end - - - - L1=t_rad(:,2:end)>t_rad(:,1); - L2=t_rad(:,1:end-1)t_rad(:,1); - [I1r,I2r] = find(L1r.*L2r); - Ir = find(L1r.*L2r); - - zir = zeros(size(Ir)); rir=zir; - indr = cell(max(I1r),1); - -for i = 1:length(Ir) - indr{I1r(i)}=[indr{I1r(i)},i]; - rir(i) = t_rad(I1r(i),1); -end - -zir = t_z(Ir) + (t_z(Ir+s1) - t_z(Ir)) ./ (t_rad(Ir+s1) - t_rad(Ir)).* (t_rad(Ir+s1) - rir) ; - - - - - av = nan(max(I1),1); - dif= nan(max(I1),1); - avs = nan(max(I1),1); - difs= nan(max(I1),1); -for i=1:length(indr) - r = zir(indr{i}); - l = zi(ind{i}); - r_i = I2r(indr{i}); - l_i = I2(ind{i}); - - m = min([length(l),length(r)]); - ix = [l(1:m),r(1:m)]'; - ix = ix(:); - - -av(i) = mean([diff(r)./diff(r_i);diff(l)./diff(l_i)]); -avs(i) = std([(r(2:end)-r(1:end-1))./(r_i(2:end)-r_i(1:end-1));(l(2:end)-l(1:end-1))./(l_i(2:end)-l_i(1:end-1))]); -dif(i) = mean((l(1:m)-r(1:m))); -difs(i) = std((l(1:m)-r(1:m))); -end - - -disp(['Larmor radius: ',num2str(mean(dif*100),'%.4f'), ' +/- ',num2str(std(difs*100),'%.4f'),' cm']) - -disp(['Av. drift: ',num2str(mean(av*1e6),'%.4f'), ' +/- ',num2str(std(avs*1e6),'%.4f'),' um/step']) - - -B = 1; -dB = 0.002 ; -dt = 1e-9; -E = 10 - -xgradB = mean(60 * 1000 ./ (B +dB*t_rad(:,1)).^ 2 *dB *dt); -xE = - E / B * dt; -disp(['expected gradB drift: ',num2str(xgradB*1e6,'%.4f'),' um/step']) -disp(['expected E field drift: ',num2str(xE*1e6,'%.4f'),' um/step']) - -if true - close all - figure - plot(t_rad(1,:),t_z(1,:)) - hold on - plot(t_rad(1,1)*ones(size(ind{1})),zi(ind{1}),'rx') -%plot(t_rad(I(indr{1})),t_z(I(indr{1})),'rx') - % plot(t_rad(1,1)*ones(size(indr{1})),zir(indr{1}),'gx') -end - diff --git a/plotter/plot3dtaiga.m b/plotter/plot3dtaiga.m deleted file mode 100644 index 3e749e5..0000000 --- a/plotter/plot3dtaiga.m +++ /dev/null @@ -1,28 +0,0 @@ -% Plot TAIGA trajectories - - -mainfolder = '../results'; -energy=60; -shotnumber = '11344_0285' -runnumber = '06May2016_140034' - - - -load([mainfolder,'/',shotnumber,'/',runnumber,'/t_rad.dat']) -load([mainfolder,'/',shotnumber,'/',runnumber,'/t_z.dat']) -load([mainfolder,'/',shotnumber,'/',runnumber,'/t_tor.dat']) - - -[R,Z]=meshgrid(rcord,zcord); - -%contour(R,Z,psi2',30) -axis equal -hold on -%plot(t_rad,t_z,'r') - -t_rad2 = sqrt(1+(t_tor./t_rad).^2).*t_rad; -t_rad -t_rad2 -plot3(t_rad,t_z,t_tor,'r') - -axis equal diff --git a/plotter/traj_plotter.m b/plotter/traj_plotter.m deleted file mode 100644 index 9be6966..0000000 --- a/plotter/traj_plotter.m +++ /dev/null @@ -1,58 +0,0 @@ -function traj_plotter(varargin) - - - mainfolder = '../results'; - - - if nargin >= 1 - shotnumber = varargin{1}; - else - shotnumber = '11774_1000'; - - end - if nargin >= 2 - runnumber = varargin{2}; - else - runnumber = '0'; - end - - - load([mainfolder,'/',shotnumber,'/',runnumber,'/t_rad.dat']); - load([mainfolder,'/',shotnumber,'/',runnumber,'/t_z.dat']); - load([mainfolder,'/',shotnumber,'/',runnumber,'/t_tor.dat']); - - - if nargin >= 3 - t_index = varargin{2}; - t_rad = t_rad(:,t_index); - t_z = t_z (:,t_index); - t_tor = t_tor(:,t_index); - end - - figure - subplot(2,2,1) - plot(t_tor,t_z) - xlabel('tor') - ylabel('z') - axis equal - - subplot(2,2,2) - plot(t_rad,t_z) - xlabel('R') - ylabel('z') - axis equal - - subplot(2,2,3) - plot(t_tor,t_rad) - xlabel('tor') - ylabel('R') - axis equal - - subplot(2,2,4) - plot3(t_rad,t_z,t_tor) - xlabel('R') - ylabel('z') - zlabel('tor') - axis equal - -end From c0e5c7811a6aa3b2ffec7cf1bc3d7d34fade3f38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Fri, 24 May 2024 00:05:45 +0200 Subject: [PATCH 12/15] Refactor plotters --- plotter/detector.py | 85 +++++++------- plotter/detector_plane.py | 170 ++++++++++++--------------- plotter/traj_plotter.py | 235 ++++++++++++++++++-------------------- 3 files changed, 229 insertions(+), 261 deletions(-) diff --git a/plotter/detector.py b/plotter/detector.py index d53c4cf..7827ca6 100644 --- a/plotter/detector.py +++ b/plotter/detector.py @@ -1,65 +1,62 @@ import sys +import os import numpy as np -import matplotlib.pyplot as pl +import matplotlib.pyplot as plt import matplotlib.patches as patches import matplotlib.cm as colormap import scipy.stats as st def detector(shotnumber, time, runnumber, is_debug=False): + results_folder = os.path.join('results', f'{shotnumber}_{time}', runnumber) - results_folder = 'results/'+shotnumber+'_'+time+'/'+runnumber+'/' - try: - C = np.loadtxt(results_folder+'detector/cellcounter.dat') - X = np.loadtxt(results_folder+'detector/detx') - Y = np.loadtxt(results_folder+'detector/dety') - except: - print ('Invalid input folder: '+results_folder) - - - xmin = X[0,0] - xmax = X[-1,-1] - ymin = Y[0,0] - ymax = Y[-1,-1] - xmean = (X[:,0]+X[:,1])/2 - ymean = (Y[:,0]+Y[:,1])/2 - cmax = C.max() + cell_counts = np.loadtxt(os.path.join(results_folder, 'detector', 'cellcounter.dat')) + det_x = np.loadtxt(os.path.join(results_folder, 'detector', 'detx')) + det_y = np.loadtxt(os.path.join(results_folder, 'detector', 'dety')) + except FileNotFoundError: + print(f'Invalid input folder: {results_folder}') + return + + xmin, xmax = det_x[0, 0], det_x[-1, -1] + ymin, ymax = det_y[0, 0], det_y[-1, -1] + xmean = (det_x[:, 0] + det_x[:, 1]) / 2 + ymean = (det_y[:, 0] + det_y[:, 1]) / 2 + cmax = cell_counts.max() if cell_counts.max() > 0 else 1 - if cmax == 0: - cmax = 1 - try: - fig = pl.figure() - ax = fig.gca() + fig, ax = plt.subplots() ax.set_xlim(xmin, xmax) ax.set_ylim(ymin, ymax) ax.set_aspect('equal') - + # Contourf plot - pl.contourf(xmean, ymean, C, cmap='afmhot', corner_mask=False) - + plt.contourf(xmean, ymean, cell_counts, cmap='afmhot', corner_mask=False) + # Detector cell matrix plot - for i in range (X[:,0].size): - for j in range (Y[:,0].size): - rect = patches.Rectangle((X[i,0], Y[j,0]), X[i,1]-X[i,0], Y[j,1]-Y[j,0], linewidth=1, - edgecolor='k', facecolor=colormap.afmhot(C[j,i]/cmax)) + for i in range(det_x[:, 0].size): + for j in range(det_y[:, 0].size): + rect = patches.Rectangle((det_x[i, 0], det_y[j, 0]), + det_x[i, 1] - det_x[i, 0], + det_y[j, 1] - det_y[j, 0], + linewidth=1, + edgecolor='k', + facecolor=colormap.afmhot(cell_counts[j, i] / cmax)) ax.add_patch(rect) if is_debug: - print(C) - - pl.xlabel(r"$x \mathrm{ [mm]}$") - pl.ylabel(r"$y \mathrm{ [mm]}$") - pl.title(r"COMPASS #"+shotnumber+" (t="+time+" s)") - - try: - print('Save plot to '+results_folder+'detpy2_'+shotnumber+'_'+time+'.pdf') - pl.savefig(results_folder+'detpy2_'+shotnumber+'_'+time+'.pdf') - pl.clf() - except: - print('Unable to save to : '+results_folder) - pl.show() - except: - print('Unable to plot: '+results_folder) + print(cell_counts) + + plt.xlabel(r"$x \mathrm{ [mm]}$") + plt.ylabel(r"$y \mathrm{ [mm]}$") + plt.title(f"COMPASS #{shotnumber} (t={time} s)") + try: + print(f'Save plot to {os.path.join(results_folder, f"detpy2_{shotnumber}_{time}.pdf")}') + plt.savefig(os.path.join(results_folder, f"detpy2_{shotnumber}_{time}.pdf")) + plt.clf() + except Exception as e: + print(f'Unable to save to: {results_folder} ({e})') + plt.show() + except Exception as e: + print(f'Unable to plot: {results_folder} ({e})') diff --git a/plotter/detector_plane.py b/plotter/detector_plane.py index 34ddf7b..99fccf0 100644 --- a/plotter/detector_plane.py +++ b/plotter/detector_plane.py @@ -1,101 +1,81 @@ import sys +import os import numpy as np -import matplotlib.pyplot as pl +import matplotlib.pyplot as plt import scipy.stats as st def detector_plane(shotnumber, time, runnumber, detector_par): - - max_distance = 1e-4 - - detector_par_list = detector_par.split(',') - - try: - R0 = float(detector_par_list[0]); - except: - R0 = 0.685 - - try: - Z0 = float(detector_par_list[1]); - except: - Z0 = 0.23 - - try: - T0 = float(detector_par_list[2]); - except: - T0 = 0.0 - - try: - detector_RZ_angle = float(detector_par_list[3]); - except: - detector_RZ_angle = 38.0 - - try: - detector_ZT_angle = float(detector_par_list[4]); - except: - detector_ZT_angle = 0.0 - - results_folder = 'results/'+shotnumber+'_'+time+'/'+runnumber+'/' - - try: - R = np.loadtxt(results_folder+'rad.dat') #, delimiter="\t" - Z = np.loadtxt(results_folder+'z.dat') - T = np.loadtxt(results_folder+'tor.dat') - except: - print('Invalid input folder: ' + results_folder) - - try: - particle_on_detector = np.abs((Z-Z0) *np.tan((detector_RZ_angle/180.0)*np.pi) + (R-R0)) < max_distance - - x = T[particle_on_detector] - y = Z[particle_on_detector] - xmin, xmax = -4e-2, 4e-2 - ymin, ymax = 18e-2, 26e-2 - - # Peform the kernel density estimate - xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] - positions = np.vstack([xx.ravel(), yy.ravel()]) - values = np.vstack([x, y]) - kernel = st.gaussian_kde(values,0.05) - f = np.reshape(kernel(positions).T, xx.shape) - - fig = pl.figure() - ax = fig.gca() - ax.set_xlim(xmin, xmax) - ax.set_ylim(ymin, ymax) - ax.set_aspect('equal') - - # Contourf plot - #cfset = ax.contourf(xx, yy, f, 500, cmap='hot') # Blues, afmhot; magma, inferno - ax.hist2d(x, y, bins=250, range = [[xmin,xmax],[ymin,ymax]], cmap='afmhot') - # Contour plot - # Label plot - - pl.xlabel(r"$T \mathrm{ [m]}$") - pl.ylabel(r"$Z \mathrm{ [m]}$") - except: - print('Unable to plot: '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') - - try: - print('Save plot to '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') - pl.savefig(results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') - pl.clf() - - except: - print('Unable to save to : '+results_folder+'detpy_'+shotnumber+'_'+time+'.pdf') - pl.show() - - xmin, xmax = -4e-2, 4e-2 - ymin, ymax = -4e-2, 4e-2 - - #x = T[particle_on_detector]-T0 - #y = (Z[particle_on_detector]-Z0) / np.sin( (detector_RZ_angle/180.0)*np.pi ) - #nbins=250 - #fig = pl.figure() - #ax = fig.gca() - #ax.set_xlim(xmin, xmax) - #ax.set_ylim(ymin, ymax) - #ax.set_aspect('equal') - #cfset = ax.hist2d(x, y, bins=nbins, range = [[xmin,xmax],[ymin,ymax]], cmap='afmhot') - #pl.colorbar(cfset[3]),ax=ax) - #pl.show() + max_distance = 1e-4 + detector_par_list = detector_par.split(',') + + try: + R0 = float(detector_par_list[0]) + except ValueError: + R0 = 0.685 + + try: + Z0 = float(detector_par_list[1]) + except ValueError: + Z0 = 0.23 + + try: + T0 = float(detector_par_list[2]) + except ValueError: + T0 = 0.0 + + try: + detector_RZ_angle = float(detector_par_list[3]) + except ValueError: + detector_RZ_angle = 38.0 + + try: + detector_ZT_angle = float(detector_par_list[4]) + except ValueError: + detector_ZT_angle = 0.0 + + results_folder = os.path.join('results', f'{shotnumber}_{time}', runnumber) + + try: + R = np.loadtxt(os.path.join(results_folder, 'rad.dat')) + Z = np.loadtxt(os.path.join(results_folder, 'z.dat')) + T = np.loadtxt(os.path.join(results_folder, 'tor.dat')) + except FileNotFoundError: + print(f'Invalid input folder: {results_folder}') + return + + try: + particle_on_detector = np.abs((Z - Z0) * np.tan((detector_RZ_angle / 180.0) * np.pi) + (R - R0)) < max_distance + x = T[particle_on_detector] + y = Z[particle_on_detector] + xmin, xmax = -4e-2, 4e-2 + ymin, ymax = 18e-2, 26e-2 + + # Perform the kernel density estimate + xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] + positions = np.vstack([xx.ravel(), yy.ravel()]) + values = np.vstack([x, y]) + kernel = st.gaussian_kde(values, 0.05) + f = np.reshape(kernel(positions).T, xx.shape) + + fig, ax = plt.subplots() + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + ax.set_aspect('equal') + + # Create a 2D histogram plot + ax.hist2d(x, y, bins=250, range=[[xmin, xmax], [ymin, ymax]], cmap='afmhot') + + # Set labels + plt.xlabel(r"$T \mathrm{ [m]}$") + plt.ylabel(r"$Z \mathrm{ [m]}$") + + try: + print(f'Save plot to {os.path.join(results_folder, f"detpy_{shotnumber}_{time}.pdf")}') + plt.savefig(os.path.join(results_folder, f"detpy_{shotnumber}_{time}.pdf")) + plt.clf() + except Exception as e: + print(f'Unable to save to: {results_folder} ({e})') + plt.show() + except Exception as e: + print(f'Unable to plot: {results_folder} ({e})') diff --git a/plotter/traj_plotter.py b/plotter/traj_plotter.py index c5a2a24..9bcd738 100644 --- a/plotter/traj_plotter.py +++ b/plotter/traj_plotter.py @@ -1,169 +1,160 @@ import sys import os -import numpy +import numpy as np from scipy.interpolate import interp1d import h5py import matplotlib.pyplot as plt - def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy): - rootFolder = os.path.abspath(os.getcwd()) - resultFolder = rootFolder+'/results' - fieldFolder = rootFolder+'/input/fieldGrid' - cdbFolder = rootFolder+'/input/cdb' - ionProfileFolder =rootFolder+'/input/ionProf' - - shotAndTime = shotnumber + '_' + time + root_folder = os.path.abspath(os.getcwd()) + result_folder = os.path.join(root_folder, 'results') + field_folder = os.path.join(root_folder, 'input', 'fieldGrid') + cdb_folder = os.path.join(root_folder, 'input', 'cdb') + ion_profile_folder = os.path.join(root_folder, 'input', 'ionProf') - workingFolder = resultFolder+'/'+shotAndTime+'/'+runnumber + shot_and_time = f"{shotnumber}_{time}" + working_folder = os.path.join(result_folder, shot_and_time, runnumber) - #ionisation - ionR = numpy.genfromtxt(fname=ionProfileFolder+'/'+shotAndTime+'/'+species+'/'+energy+'/rad.dat') - ionY = numpy.genfromtxt(fname=ionProfileFolder+'/'+shotAndTime+'/'+species+'/'+energy+'/ionyeald.dat') - ionX1 = numpy.gradient(ionY) - ionF = interp1d(ionR, -ionX1/numpy.amax(numpy.absolute(ionX1))) + # Ionization + ionR = np.genfromtxt(fname=os.path.join(ion_profile_folder, shot_and_time, species, energy, 'rad.dat')) + ionY = np.genfromtxt(fname=os.path.join(ion_profile_folder, shot_and_time, species, energy, 'ionyeald.dat')) + ionX1 = np.gradient(ionY) + ionF = interp1d(ionR, -ionX1 / np.amax(np.absolute(ionX1))) plt.plot(ionR, ionF(ionR)) try: - print('Save plot to '+workingFolder+'/ion_'+shotAndTime+'.pdf') - plt.savefig(workingFolder+'/ion_'+shotAndTime+'.pdf') + print(f"Save plot to {os.path.join(working_folder, 'ion_' + shot_and_time + '.pdf')}") + plt.savefig(os.path.join(working_folder, f"ion_{shot_and_time}.pdf")) plt.clf() - - except: - print('Unable to save to : '+workingFolder+'/ion_'+shotAndTime+'.pdf') + except Exception as e: + print(f"Unable to save to: {working_folder} ({e})") plt.show() - # trajectories - t_rad = numpy.genfromtxt(fname=workingFolder+'/t_rad.dat') - t_z = numpy.genfromtxt(fname=workingFolder+'/t_z.dat') - t_tor = numpy.genfromtxt(fname=workingFolder+'/t_tor.dat') - cellid = numpy.genfromtxt(fname=workingFolder+'/detector/cellid.dat') - - rcord = numpy.genfromtxt(fname=fieldFolder+'/'+shotAndTime+'/rcord.dat') - zcord = numpy.genfromtxt(fname=fieldFolder+'/'+shotAndTime+'/zcord.dat') - psi2 = numpy.genfromtxt(fname=fieldFolder+'/'+shotAndTime+'/psi2.dat') - - hdf5Path = cdbFolder+'/'+shotnumber+'/EFITXX/EFITXX.1.h5' - hdf5File = h5py.File(hdf5Path, 'r') - hdf5Times = hdf5File.get('time')[()] - hdf5PsiBoundary = hdf5File.get('output/globalParameters/psiBoundary')[()] - hdf5PsiAxis = hdf5File.get('output/globalParameters/psiAxis')[()] - hdf5BoundaryR = hdf5File.get('output/separatrixGeometry/boundaryCoordsR')[()] - hdf5BoundaryZ = hdf5File.get('output/separatrixGeometry/boundaryCoordsZ')[()] - - hdf5Index = int(numpy.argwhere(hdf5Times == float(time)/1000.)) - - psiAxis = hdf5PsiAxis[hdf5Index] - psiBoundary = hdf5PsiBoundary[hdf5Index] - boundaryR = hdf5BoundaryR[hdf5Index,:] - boundaryZ = hdf5BoundaryZ[hdf5Index,:] - - PSI = numpy.transpose(psi2)/psiAxis - - R, Z = numpy.meshgrid(rcord, zcord) - + # Trajectories + t_rad = np.genfromtxt(fname=os.path.join(working_folder, 't_rad.dat')) + t_z = np.genfromtxt(fname=os.path.join(working_folder, 't_z.dat')) + t_tor = np.genfromtxt(fname=os.path.join(working_folder, 't_tor.dat')) + cellid = np.genfromtxt(fname=os.path.join(working_folder, 'detector', 'cellid.dat')) + + rcord = np.genfromtxt(fname=os.path.join(field_folder, shot_and_time, 'rcord.dat')) + zcord = np.genfromtxt(fname=os.path.join(field_folder, shot_and_time, 'zcord.dat')) + psi2 = np.genfromtxt(fname=os.path.join(field_folder, shot_and_time, 'psi2.dat')) + + hdf5_path = os.path.join(cdb_folder, shotnumber, 'EFITXX', 'EFITXX.1.h5') + hdf5_file = h5py.File(hdf5_path, 'r') + hdf5_times = hdf5_file.get('time')[()] + hdf5_psi_boundary = hdf5_file.get('output/globalParameters/psiBoundary')[()] + hdf5_psi_axis = hdf5_file.get('output/globalParameters/psiAxis')[()] + hdf5_boundary_r = hdf5_file.get('output/separatrixGeometry/boundaryCoordsR')[()] + hdf5_boundary_z = hdf5_file.get('output/separatrixGeometry/boundaryCoordsZ')[()] + hdf5_index = int(np.argwhere(hdf5_times == float(time) / 1000.)) + psi_axis = hdf5_psi_axis[hdf5_index] + psi_boundary = hdf5_psi_boundary[hdf5_index] + boundary_r = hdf5_boundary_r[hdf5_index, :] + boundary_z = hdf5_boundary_z[hdf5_index, :] + + # Calculate PSI + PSI = np.transpose(psi2) / psi_axis + + # Create meshgrid + R, Z = np.meshgrid(rcord, zcord) + + # Filter data t_index = cellid >= 0 - t_rad = t_rad[:,t_index] - t_z = t_z [:,t_index] - t_tor = t_tor[:,t_index] - numberOfIons = numpy.size(t_rad, 1) - - #ion histogram (real) - fig = plt.figure() - plt.hist(t_rad[0,:], density=True, bins=50) - plt.xlabel(r'$R$ [m]') - plt.title (r'COMPASS #'+shotnumber+'\n ($t = $'+time+' ms)') - + t_rad = t_rad[:, t_index] + t_z = t_z[:, t_index] + t_tor = t_tor[:, t_index] + number_of_ions = np.size(t_rad, 1) + + # Plot ion histogram + fig, ax = plt.subplots() + ax.hist(t_rad[0, :], density=True, bins=50) + ax.set_xlabel(r'$R$ [m]') + ax.set_title(f'COMPASS #{shotnumber}\n($t = {time}$ ms)') + + # Save plot + plot_filename = os.path.join(working_folder, f'start_{shot_and_time}.pdf') try: - print('Save plot to '+workingFolder+'/start_'+shotAndTime+'.pdf') - plt.savefig(workingFolder+'/start_'+shotAndTime+'.pdf') + plt.savefig(plot_filename) plt.clf() - + print(f'Saved plot to {plot_filename}') except: - print('Unable to save to : '+workingFolder+'/start_'+shotAndTime+'.pdf') + print(f'Unable to save to: {plot_filename}') plt.show() - #detector - plt.clf() + # Create detector plots fig, axs = plt.subplots(2, 2) - axs[0,0].plot(t_rad[-1,:], t_z[-1,:], '.') - #axs[0,0].set_xlabel(r'$R$ [m]') - axs[0,0].set_ylabel(r'$Z$ [m]') - axs[0,0].yaxis.set_ticks_position('both') + axs[0, 0].plot(t_rad[-1, :], t_z[-1, :], '.') + axs[0, 0].set_ylabel(r'$Z$ [m]') + axs[0, 0].yaxis.set_ticks_position('both') - axs[0,1].plot(t_tor[-1,:], t_z[-1,:], '.') - axs[0,1].set_xlabel(r'$T$ [m]') - axs[0,1].set_ylabel(r'$Z$ [m]') - axs[0,1].yaxis.tick_right() - axs[0,1].yaxis.set_ticks_position('both') + axs[0, 1].plot(t_tor[-1, :], t_z[-1, :], '.') + axs[0, 1].set_xlabel(r'$T$ [m]') + axs[0, 1].set_ylabel(r'$Z$ [m]') + axs[0, 1].yaxis.tick_right() + axs[0, 1].yaxis.set_ticks_position('both') - axs[1,0].plot(t_rad[-1,:], t_tor[-1,:], '.') - axs[1,0].set_xlabel(r'$R$ [m]') - axs[1,0].set_ylabel(r'$T$ [m]') - axs[1,0].yaxis.set_ticks_position('both') + axs[1, 0].plot(t_rad[-1, :], t_tor[-1, :], '.') + axs[1, 0].set_xlabel(r'$R$ [m]') + axs[1, 0].set_ylabel(r'$T$ [m]') + axs[1, 0].yaxis.set_ticks_position('both') - axs[1,1].axis('off') + axs[1, 1].axis('off') - fig.suptitle(r'COMPASS #'+shotnumber+' ($t = $'+time+' ms)') - #plt.colorbar() + fig.suptitle(f'COMPASS #{shotnumber} ($t = {time}$ ms)') + plt.show() - try: - print('Save plot to '+workingFolder+'/end_'+shotAndTime+'.pdf') - plt.savefig(workingFolder+'/end_'+shotAndTime+'.pdf') - plt.clf() - - except: - print('Unable to save to : '+workingFolder+'/end_'+shotAndTime+'.pdf') - plt.show() - - #2D figure + # Create a 2D figure plt.clf() fig = plt.figure() ax = fig.add_subplot(111) - #plt.contourf(R, Z, -PSI, levels=numpy.arange(-1.05,2,.05), cmap='Spectral') - #plt.contourf(R, Z, PSI, levels=numpy.arange(-2,2,.05), cmap='coolwarm') - #plt.contourf(R, Z, PSI, levels=numpy.arange(-2,2,.05), cmap='twilight_shifted') - plt.contourf(R, Z, PSI, levels=numpy.arange(-2,2,.1), cmap='GnBu') - plt.plot(boundaryR, boundaryZ, color=(0,0,0)) - for i in range(numberOfIons): - i1 = 1 - i2 = numpy.asscalar(ionF(t_rad[1,i])) - plt.plot([rcord[-1], t_rad[1,i]], [t_z[0,i], t_z[0,i]], color=(i1, i1, 0)) - if i2<0: - i2=0 - if i2>1: - i2=1 - plt.plot(t_rad[:,i], t_z[:,i], color=(i2, 0, 0)) - - # plot detector + + # Plot filled contours + plt.contourf(R, Z, PSI, levels=np.arange(-2, 2, 0.1), cmap='GnBu') + + # Plot boundary line + plt.plot(boundary_r, boundary_z, color='black') + + # Iterate over ions + for i in range(number_of_ions): + ion_intensity = 1 + ion_factor = np.asscalar(ionF(t_rad[1, i])) + + # Adjust ion factor within bounds + ion_factor = max(0, min(1, ion_factor)) + + # Plot ion trajectory + plt.plot([rcord[-1], t_rad[1, i]], [t_z[0, i], t_z[0, i]], color=(ion_intensity, ion_intensity, 0)) + plt.plot(t_rad[:, i], t_z[:, i], color=(ion_factor, 0, 0)) + try: + # Detector parameters detector_y_size = 12.2e-3 detector_x0, detector_y0, detector_z0, detector_phi, detector_psi = map(float, detector_par.split(',')) - detector_x_tilt = detector_y_size * numpy.sin(detector_phi / 180 * numpy.pi) - detector_y_tilt = detector_y_size * numpy.cos(detector_phi / 180 * numpy.pi) - plt.plot([detector_x0, detector_x0], [detector_y0, max(zcord)], color=(0, 0, 1)) - plt.plot([detector_x0, detector_x0], [detector_y0, max(zcord)], color=(0, 0, 1)) + # Calculate tilted coordinates + detector_x_tilt = detector_y_size * np.sin(detector_phi / 180 * np.pi) + detector_y_tilt = detector_y_size * np.cos(detector_phi / 180 * np.pi) + + # Plot detector lines + plt.plot([detector_x0, detector_x0], [detector_y0, max(zcord)], + color='blue', linewidth=1) plt.plot([detector_x0 - detector_x_tilt, detector_x0 + detector_x_tilt], [detector_y0 + detector_y_tilt, detector_y0 - detector_y_tilt], - color=(0, 0, 1), linewidth=2) - except: + color='blue', linewidth=2) + except ValueError: print('Invalid detector coordinates.') - pass plt.xlabel(r'$R$ [m]') plt.ylabel(r'$Z$ [m]') plt.title (r'COMPASS #'+shotnumber+'\n ($t = $'+time+' ms)') ax.set_aspect('equal', 'box') - #plt.colorbar() - try: - print('Save plot to '+workingFolder+'/traj_'+shotAndTime+'.pdf') - plt.savefig(workingFolder+'/traj_'+shotAndTime+'.pdf') + plot_filename = os.path.join(working_folder, f'traj_{shot_and_time}.pdf') + plt.savefig(plot_filename) plt.clf() - + print(f'Saved plot to {plot_filename}') except: - print('Unable to save to : '+workingFolder+'/traj_'+shotAndTime+'.pdf') + print(f'Unable to save to: {plot_filename}') plt.show() - From bf5b59e3cb8a2182ae5aad16589e7242b295cb20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Fri, 24 May 2024 00:07:17 +0200 Subject: [PATCH 13/15] Refactor plotters --- plotter/traj_plotter.py | 1 + 1 file changed, 1 insertion(+) diff --git a/plotter/traj_plotter.py b/plotter/traj_plotter.py index 9bcd738..2dc3bc0 100644 --- a/plotter/traj_plotter.py +++ b/plotter/traj_plotter.py @@ -5,6 +5,7 @@ import h5py import matplotlib.pyplot as plt + def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy): root_folder = os.path.abspath(os.getcwd()) result_folder = os.path.join(root_folder, 'results') From 47329fb0e72268dbde22d7821ea6fded35817874 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Wed, 29 May 2024 23:21:07 +0200 Subject: [PATCH 14/15] Refactor plotters --- plotter/traj_plotter.py | 43 +++++++++++++++++++++++++---------------- sd.py | 13 ++++++++----- 2 files changed, 34 insertions(+), 22 deletions(-) diff --git a/plotter/traj_plotter.py b/plotter/traj_plotter.py index 2dc3bc0..3cb3f89 100644 --- a/plotter/traj_plotter.py +++ b/plotter/traj_plotter.py @@ -1,12 +1,14 @@ import sys import os import numpy as np +import numpy.matlib from scipy.interpolate import interp1d import h5py import matplotlib.pyplot as plt -def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy): + +def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, title=None): root_folder = os.path.abspath(os.getcwd()) result_folder = os.path.join(root_folder, 'results') field_folder = os.path.join(root_folder, 'input', 'fieldGrid') @@ -16,6 +18,9 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy): shot_and_time = f"{shotnumber}_{time}" working_folder = os.path.join(result_folder, shot_and_time, runnumber) + if title is None: + r'COMPASS #' + shotnumber + '\n ($t = $' + time + ' ms)' + # Ionization ionR = np.genfromtxt(fname=os.path.join(ion_profile_folder, shot_and_time, species, energy, 'rad.dat')) ionY = np.genfromtxt(fname=os.path.join(ion_profile_folder, shot_and_time, species, energy, 'ionyeald.dat')) @@ -72,14 +77,14 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy): fig, ax = plt.subplots() ax.hist(t_rad[0, :], density=True, bins=50) ax.set_xlabel(r'$R$ [m]') - ax.set_title(f'COMPASS #{shotnumber}\n($t = {time}$ ms)') + ax.set_title(title) # Save plot plot_filename = os.path.join(working_folder, f'start_{shot_and_time}.pdf') try: plt.savefig(plot_filename) plt.clf() - print(f'Saved plot to {plot_filename}') + print(f'Save plot to {plot_filename}') except: print(f'Unable to save to: {plot_filename}') plt.show() @@ -103,7 +108,7 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy): axs[1, 1].axis('off') - fig.suptitle(f'COMPASS #{shotnumber} ($t = {time}$ ms)') + fig.suptitle(title) plt.show() # Create a 2D figure @@ -118,17 +123,19 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy): plt.plot(boundary_r, boundary_z, color='black') # Iterate over ions - for i in range(number_of_ions): - ion_intensity = 1 - ion_factor = np.asscalar(ionF(t_rad[1, i])) - + intensity_indices = np.argsort(ionF(t_rad[1, :])) + ion_factor = ionF(t_rad[1, :]) + for i in intensity_indices: # Adjust ion factor within bounds - ion_factor = max(0, min(1, ion_factor)) - + ion_factor[i] = max(0, min(1, ion_factor[i])) # Plot ion trajectory - plt.plot([rcord[-1], t_rad[1, i]], [t_z[0, i], t_z[0, i]], color=(ion_intensity, ion_intensity, 0)) - plt.plot(t_rad[:, i], t_z[:, i], color=(ion_factor, 0, 0)) + plt.plot(t_rad[:, i], t_z[:, i], color=(ion_factor[i], 0, 0)) + alpha_factor = 50/number_of_ions + for i in intensity_indices: + # Plot atom trajectory + plt.plot([rcord[-1], t_rad[1, i]], [t_z[0, i], t_z[0, i]], color=(1, 1, 0), alpha=ion_factor[i] * alpha_factor) + print(number_of_ions) try: # Detector parameters detector_y_size = 12.2e-3 @@ -149,13 +156,15 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy): plt.xlabel(r'$R$ [m]') plt.ylabel(r'$Z$ [m]') - plt.title (r'COMPASS #'+shotnumber+'\n ($t = $'+time+' ms)') + plt.title(title) ax.set_aspect('equal', 'box') try: - plot_filename = os.path.join(working_folder, f'traj_{shot_and_time}.pdf') - plt.savefig(plot_filename) + plot_filename = os.path.join(working_folder, f'traj_{shot_and_time}') + plt.savefig(plot_filename+'.pdf', dpi=300, bbox_inches='tight') + print(f'Save plot to {plot_filename}.pdf') + plt.savefig(plot_filename+'.png', dpi=300, bbox_inches='tight') + print(f'Save plot to {plot_filename}.png') plt.clf() - print(f'Saved plot to {plot_filename}') except: - print(f'Unable to save to: {plot_filename}') + print(f'Unable to save to: {plot_filename}.pdf/png') plt.show() diff --git a/sd.py b/sd.py index 53f1e05..2aed952 100644 --- a/sd.py +++ b/sd.py @@ -11,7 +11,7 @@ class SD: def __init__(self, shot_number, time, species, energy, detector_par, - runnumber=None, is_trajectory_detailed=False): + runnumber=None, is_trajectory_detailed=False, title=None): self.shot_number = str(int(shot_number)) self.time = str(int(time)) self.species = species @@ -19,6 +19,7 @@ def __init__(self, shot_number, time, species, energy, detector_par, self.detector_par = detector_par self.runnumber = runnumber self.is_trajectory_detailed = is_trajectory_detailed + self.title = title self.is_run_simulation = True self.set_runnumber() if self.is_run_simulation: @@ -48,7 +49,7 @@ def run_plotters(self): detector(self.shot_number, self.time, self.runnumber) detector_plane(self.shot_number, self.time, self.runnumber, self.detector_par) if self.is_trajectory_detailed: - traj_plotter(self.shot_number, self.time, self.runnumber, self.detector_par, self.species, self.energy) + traj_plotter(self.shot_number, self.time, self.runnumber, self.detector_par, self.species, self.energy, self.title) def set_runnumber(self): if self.runnumber is None: @@ -65,7 +66,7 @@ def write_parameter_file(self): f.write("energy=" + self.energy + " #keV\n") f.write("toroidal_deflection=0 #degree; + 200 V deflection\n") f.write("diameter=5 #mm\n") - f.write("particles=1000\n") + f.write("particles=10000\n") f.write("detector='" + self.detector_par + "'\n") f.write("electric_field_module=0\n") f.write("detector_mask='final'\n") @@ -82,7 +83,9 @@ def delete_parameter_file(self): a_time = 1097 a_species = 'Li' an_energy = 70 - a_detector = "0.6846,0.253,0.0,38,0" + z_det = 0.253 + a_detector = f"0.6846,{z_det},0.0,38,0" a_runnumber = '20240523214835' is_plot_trajectory = True - SD(a_shot_number, a_time, a_species, an_energy, a_detector, a_runnumber, is_trajectory_detailed=is_plot_trajectory) + a_title = f"Reference discharge\n{a_species} beam, $E ={an_energy}$ keV\n$Z_D={z_det*100}$ cm" + SD(a_shot_number, a_time, a_species, an_energy, a_detector, a_runnumber, is_trajectory_detailed=is_plot_trajectory, title=a_title) From 4e7426ea013da57f91eb4e75d45c449cbb522d8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?M=C3=A1ty=C3=A1s=20Aradi?= <20679946+matyasaradi@users.noreply.github.com> Date: Fri, 2 Aug 2024 02:17:17 +0200 Subject: [PATCH 15/15] Update --- plotter/detector.py | 26 +++++++++--- plotter/detector_plane.py | 25 ++++++++++-- plotter/traj_plotter.py | 25 ++++++++---- preproc/renate_od/shake_flux_surface.py | 53 +++++++++++++++++-------- sd.py | 15 +++---- 5 files changed, 102 insertions(+), 42 deletions(-) diff --git a/plotter/detector.py b/plotter/detector.py index 7827ca6..7f697f3 100644 --- a/plotter/detector.py +++ b/plotter/detector.py @@ -5,11 +5,15 @@ import matplotlib.patches as patches import matplotlib.cm as colormap import scipy.stats as st +import qrcode -def detector(shotnumber, time, runnumber, is_debug=False): +def detector(shotnumber, time, runnumber, title=None, is_debug=False): results_folder = os.path.join('results', f'{shotnumber}_{time}', runnumber) + if title is None: + r'COMPASS #' + shotnumber + '\n ($t = $' + time + ' ms)' + try: cell_counts = np.loadtxt(os.path.join(results_folder, 'detector', 'cellcounter.dat')) det_x = np.loadtxt(os.path.join(results_folder, 'detector', 'detx')) @@ -25,6 +29,8 @@ def detector(shotnumber, time, runnumber, is_debug=False): cmax = cell_counts.max() if cell_counts.max() > 0 else 1 try: + plt.rc('font', family='serif') + plt.rc('text', usetex=True) fig, ax = plt.subplots() ax.set_xlim(xmin, xmax) ax.set_ylim(ymin, ymax) @@ -47,13 +53,21 @@ def detector(shotnumber, time, runnumber, is_debug=False): if is_debug: print(cell_counts) - plt.xlabel(r"$x \mathrm{ [mm]}$") - plt.ylabel(r"$y \mathrm{ [mm]}$") - plt.title(f"COMPASS #{shotnumber} (t={time} s)") + plt.xlabel(r"toroidal detector position [mm]") + plt.ylabel(r"vertical detector position [mm]") + plt.title(title) + + qr_content = f'raw.githubusercontent.com/taiga-project/refs/run/{runnumber}' + print(f'Add QR code: {qr_content}') + qr_image = qrcode.make(qr_content) + ax_qr = fig.add_axes([0.175, 0.495, 0.08, 0.5], anchor='NE', zorder=10) + ax_qr.imshow(qr_image, cmap='gray') + ax_qr.axis('off') try: - print(f'Save plot to {os.path.join(results_folder, f"detpy2_{shotnumber}_{time}.pdf")}') - plt.savefig(os.path.join(results_folder, f"detpy2_{shotnumber}_{time}.pdf")) + filename = f"detector_{runnumber}.pdf" + print(f'Save plot to {os.path.join(results_folder, filename)}') + plt.savefig(os.path.join(results_folder, filename), dpi=300, bbox_inches='tight') plt.clf() except Exception as e: print(f'Unable to save to: {results_folder} ({e})') diff --git a/plotter/detector_plane.py b/plotter/detector_plane.py index 99fccf0..c4fa377 100644 --- a/plotter/detector_plane.py +++ b/plotter/detector_plane.py @@ -3,9 +3,10 @@ import numpy as np import matplotlib.pyplot as plt import scipy.stats as st +import qrcode -def detector_plane(shotnumber, time, runnumber, detector_par): +def detector_plane(shotnumber, time, runnumber, detector_par, title=None): max_distance = 1e-4 detector_par_list = detector_par.split(',') @@ -36,6 +37,9 @@ def detector_plane(shotnumber, time, runnumber, detector_par): results_folder = os.path.join('results', f'{shotnumber}_{time}', runnumber) + if title is None: + r'COMPASS #' + shotnumber + '\n ($t = $' + time + ' ms)' + try: R = np.loadtxt(os.path.join(results_folder, 'rad.dat')) Z = np.loadtxt(os.path.join(results_folder, 'z.dat')) @@ -50,6 +54,9 @@ def detector_plane(shotnumber, time, runnumber, detector_par): y = Z[particle_on_detector] xmin, xmax = -4e-2, 4e-2 ymin, ymax = 18e-2, 26e-2 + #ymin, ymax = 24e-2, 32e-2 + #xmin, xmax = -8e-2, 0e-2 + #ymin, ymax = 14e-2, 22e-2 # Perform the kernel density estimate xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] @@ -59,6 +66,8 @@ def detector_plane(shotnumber, time, runnumber, detector_par): f = np.reshape(kernel(positions).T, xx.shape) fig, ax = plt.subplots() + plt.rc('font', family='serif') + plt.rc('text', usetex=True) ax.set_xlim(xmin, xmax) ax.set_ylim(ymin, ymax) ax.set_aspect('equal') @@ -67,12 +76,20 @@ def detector_plane(shotnumber, time, runnumber, detector_par): ax.hist2d(x, y, bins=250, range=[[xmin, xmax], [ymin, ymax]], cmap='afmhot') # Set labels - plt.xlabel(r"$T \mathrm{ [m]}$") - plt.ylabel(r"$Z \mathrm{ [m]}$") + plt.xlabel(r"$T$ [m]") + plt.ylabel(r"$Z$ [m]") + plt.title(title) + + qr_content = f'raw.githubusercontent.com/taiga-project/refs/run/{runnumber}' + print(f'Add QR code: {qr_content}') + qr_image = qrcode.make(qr_content) + ax_qr = fig.add_axes([0.24, 0.495, 0.08, 0.5], anchor='NE', zorder=10) + ax_qr.imshow(qr_image, cmap='gray') + ax_qr.axis('off') try: print(f'Save plot to {os.path.join(results_folder, f"detpy_{shotnumber}_{time}.pdf")}') - plt.savefig(os.path.join(results_folder, f"detpy_{shotnumber}_{time}.pdf")) + plt.savefig(os.path.join(results_folder, f"detpy_{shotnumber}_{time}.pdf"), dpi=300, bbox_inches='tight') plt.clf() except Exception as e: print(f'Unable to save to: {results_folder} ({e})') diff --git a/plotter/traj_plotter.py b/plotter/traj_plotter.py index 3cb3f89..f8d9b16 100644 --- a/plotter/traj_plotter.py +++ b/plotter/traj_plotter.py @@ -5,10 +5,10 @@ from scipy.interpolate import interp1d import h5py import matplotlib.pyplot as plt +import qrcode - -def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, title=None): +def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, title=None, R_min=0.65): root_folder = os.path.abspath(os.getcwd()) result_folder = os.path.join(root_folder, 'results') field_folder = os.path.join(root_folder, 'input', 'fieldGrid') @@ -25,7 +25,7 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit ionR = np.genfromtxt(fname=os.path.join(ion_profile_folder, shot_and_time, species, energy, 'rad.dat')) ionY = np.genfromtxt(fname=os.path.join(ion_profile_folder, shot_and_time, species, energy, 'ionyeald.dat')) ionX1 = np.gradient(ionY) - ionF = interp1d(ionR, -ionX1 / np.amax(np.absolute(ionX1))) + ionF = interp1d(ionR, -ionX1 / np.amax(np.absolute(ionX1[ionR > R_min]))) plt.plot(ionR, ionF(ionR)) @@ -67,7 +67,7 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit R, Z = np.meshgrid(rcord, zcord) # Filter data - t_index = cellid >= 0 + t_index = (cellid >= 0) & (t_rad[1, :] > R_min) t_rad = t_rad[:, t_index] t_z = t_z[:, t_index] t_tor = t_tor[:, t_index] @@ -113,11 +113,14 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit # Create a 2D figure plt.clf() + plt.rc('font', family='serif') + plt.rc('text', usetex=True) fig = plt.figure() ax = fig.add_subplot(111) # Plot filled contours plt.contourf(R, Z, PSI, levels=np.arange(-2, 2, 0.1), cmap='GnBu') + #plt.contourf(R, Z, PSI, levels=np.arange(0, 12, 0.3), cmap='GnBu') # Plot boundary line plt.plot(boundary_r, boundary_z, color='black') @@ -135,7 +138,6 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit for i in intensity_indices: # Plot atom trajectory plt.plot([rcord[-1], t_rad[1, i]], [t_z[0, i], t_z[0, i]], color=(1, 1, 0), alpha=ion_factor[i] * alpha_factor) - print(number_of_ions) try: # Detector parameters detector_y_size = 12.2e-3 @@ -146,11 +148,12 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit detector_y_tilt = detector_y_size * np.cos(detector_phi / 180 * np.pi) # Plot detector lines + detector_colour = (0, 0.5, 0) plt.plot([detector_x0, detector_x0], [detector_y0, max(zcord)], - color='blue', linewidth=1) + color=detector_colour, linewidth=1) plt.plot([detector_x0 - detector_x_tilt, detector_x0 + detector_x_tilt], [detector_y0 + detector_y_tilt, detector_y0 - detector_y_tilt], - color='blue', linewidth=2) + color=detector_colour, linewidth=2) except ValueError: print('Invalid detector coordinates.') @@ -158,8 +161,14 @@ def traj_plotter(shotnumber, time, runnumber, detector_par, species, energy, tit plt.ylabel(r'$Z$ [m]') plt.title(title) ax.set_aspect('equal', 'box') + qr_content = f'https://taiga-project.github.io?{runnumber}' + print(f'Add QR code: {qr_content}') + qr_image = qrcode.make(qr_content) + ax_qr = fig.add_axes([0.24, 0.495, 0.08, 0.5], anchor='NE', zorder=10) + ax_qr.imshow(qr_image, cmap='gray') + ax_qr.axis('off') try: - plot_filename = os.path.join(working_folder, f'traj_{shot_and_time}') + plot_filename = os.path.join(working_folder, f'traj_{runnumber}') plt.savefig(plot_filename+'.pdf', dpi=300, bbox_inches='tight') print(f'Save plot to {plot_filename}.pdf') plt.savefig(plot_filename+'.png', dpi=300, bbox_inches='tight') diff --git a/preproc/renate_od/shake_flux_surface.py b/preproc/renate_od/shake_flux_surface.py index 88ac506..3071064 100644 --- a/preproc/renate_od/shake_flux_surface.py +++ b/preproc/renate_od/shake_flux_surface.py @@ -1,5 +1,8 @@ import matplotlib.pyplot import pandas +import sys +import os +sys.path.append(os.path.join(os.path.dirname(__file__), '..', '..')) from preproc.renate_od.interface import get_lcfs_radial_coordinate from preproc.renate_od.manager import RenateODManager @@ -47,7 +50,6 @@ def __init__(self, beamlet_geometry, shot_number, time, species, energy, shift=0 def shake_flux_surface(species, shot_number, time, energy): - export_directory='.' z = 0 tor = 0 @@ -58,26 +60,31 @@ def shake_flux_surface(species, shot_number, time, energy): matplotlib.pyplot.minorticks_on() matplotlib.pyplot.grid(which='major') matplotlib.pyplot.xlabel('$R$ [m]', labelpad=-10.5, loc='right') - matplotlib.pyplot.ylabel('neutral beam attenuation') - matplotlib.pyplot.title('COMPASS #' + shot_number + ' (' + time + ' ms)') - R_LCFS = get_lcfs_radial_coordinate(shot_number, time) + matplotlib.pyplot.ylabel('normalized\nneutral beam density') + matplotlib.pyplot.title('reference discharge, ' + species + ' beam, ' + str(energy) + ' keV') + R_LCFS = get_lcfs_radial_coordinate(shot_number=shot_number, time=time) matplotlib.pyplot.axvline(R_LCFS, c='red', ls='--') + matplotlib.pyplot.text(R_LCFS+0.005, 0.56, 'reference', c='red', fontsize=6) matplotlib.pyplot.text(R_LCFS+0.005, 0.45, 'LCFS', c='red', fontsize=12) export_directory = get_home_directory() + '/input/ionProf/' + shot_number + '_' + time for shift in numpy.linspace(-0.02, 0.02, 5): r = RenateODManagerMock(beamlet_geometry, shot_number, time, species, energy, shift) radial_coordinate, relative_attenuation = r.get_attenuation_profile() - ax.plot(radial_coordinate, relative_attenuation, '-', linewidth=2, label=r'$\rho$->$\rho$+%.2f'%shift) + label = r'with %3.f%% $\rho$' % (100*(1+shift)) + if shift == 0: + label = r'reference $\rho$' + ax.plot(radial_coordinate, relative_attenuation, '-', linewidth=2, label=label, + color=(max(abs(shift*10), shift*50), abs(shift*30), max(abs(shift*10), -shift*50))) relative_attenuation.fillna(0).to_csv(export_directory + '/ionyeald' + str(int(100*(1+shift)))+'.dat', index=False, header=False) - ax.legend() - matplotlib.pyplot.savefig(export_directory+'/attenuation.pdf') - matplotlib.pyplot.savefig(export_directory+'/attenuation.svg') + ax.legend(labelspacing=0.3, borderpad=0) + matplotlib.pyplot.xlim(0.6, 0.75) + matplotlib.pyplot.subplots_adjust(left=0.18, right=0.95, top=0.85, bottom=0.15) + matplotlib.pyplot.savefig(export_directory+'/attenuation_shaked.pdf') + matplotlib.pyplot.savefig(export_directory+'/attenuation_shaked.svg') def shake_flux_surface_silent(species, shot_number, time, energy): - - beamlet_geometry = BeamletGeometry() beamlet_geometry.rad = numpy.linspace(0.72, 0.6, 1000) beamlet_geometry.set_with_value(0, 'z', 'rad') @@ -95,7 +102,7 @@ def shake_flux_surface_silent(species, shot_number, time, energy): return numpy.max(diff) -if __name__ == "__main__": +def test_multi_energy(): a_species = 'Li' a_shot_number = '17178' a_time = '1097' @@ -118,9 +125,21 @@ def shake_flux_surface_silent(species, shot_number, time, energy): matplotlib.pyplot.savefig(export_directory+'/maxdiff.svg') -# if __name__ == "__main__": -# a_species = 'Li' -# a_shot_number = '17178' -# a_time = '1097' -# an_energy = 80 -# shake_flux_surface(a_species, a_shot_number, a_time, an_energy) +def test_single_energy(): + a_species = 'Li' + a_shot_number = '17178' + a_time = '1097' + an_energy = 80 + shake_flux_surface(a_species, a_shot_number, a_time, an_energy) + + +def test_single_energy2(): + a_species = 'Na' + a_shot_number = '17178' + a_time = '1097' + an_energy = 50 + shake_flux_surface(a_species, a_shot_number, a_time, an_energy) + + +if __name__ == "__main__": + test_single_energy2() diff --git a/sd.py b/sd.py index 2aed952..032b0c3 100644 --- a/sd.py +++ b/sd.py @@ -46,8 +46,8 @@ def run_core(self): def run_plotters(self): os.chdir(get_home_directory()) - detector(self.shot_number, self.time, self.runnumber) - detector_plane(self.shot_number, self.time, self.runnumber, self.detector_par) + detector(self.shot_number, self.time, self.runnumber, self.title) + detector_plane(self.shot_number, self.time, self.runnumber, self.detector_par, self.title) if self.is_trajectory_detailed: traj_plotter(self.shot_number, self.time, self.runnumber, self.detector_par, self.species, self.energy, self.title) @@ -66,7 +66,7 @@ def write_parameter_file(self): f.write("energy=" + self.energy + " #keV\n") f.write("toroidal_deflection=0 #degree; + 200 V deflection\n") f.write("diameter=5 #mm\n") - f.write("particles=10000\n") + f.write("particles=1000\n") f.write("detector='" + self.detector_par + "'\n") f.write("electric_field_module=0\n") f.write("detector_mask='final'\n") @@ -83,9 +83,10 @@ def delete_parameter_file(self): a_time = 1097 a_species = 'Li' an_energy = 70 - z_det = 0.253 - a_detector = f"0.6846,{z_det},0.0,38,0" - a_runnumber = '20240523214835' + z_det = 22 + t_det = -3.0 + a_detector = f"0.6846,{z_det/100},{t_det/100},38,0" + a_runnumber = None#'20240530015119' is_plot_trajectory = True - a_title = f"Reference discharge\n{a_species} beam, $E ={an_energy}$ keV\n$Z_D={z_det*100}$ cm" + a_title = f"reference magnetic field\n{a_species} beam, $E ={an_energy}$ keV\n$Z_D={z_det}$ cm, $T_D={t_det}$ cm" SD(a_shot_number, a_time, a_species, an_energy, a_detector, a_runnumber, is_trajectory_detailed=is_plot_trajectory, title=a_title)