polyakgergo/putdata.py at main · kdbsb1/polyakgergo · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
import numpy
from crm_solver import solve
from utility import get_data_from_hdf5
from scipy.interpolate import interpn


class Inputs:
    def __init__(self, filename='/marconi_work/eufus_gw/work/g2bszond/renate-od.git/trunk/'
                                'data/Profiles/SOL_profiles/EAST_412.h5', time_index=0):
        self.filename = filename
        self.time_index = time_index
        self.initial_condition = numpy.array([1, 0, 0, 0, 0, 0, 0, 0, 0])
        self.number_of_levels = len(self.initial_condition)
        self.step_interval = 0.08
        self.step_number = 100
        self.beam_energy = 60
        self.beam_species = 'Li'
        self.profiles = Profiles(filename, time_index=self.time_index)
        self.profiles.calculate_new_r_axis() # Shift coordinate

        self.vertical = numpy.zeros(self.step_number)
        self.steps = numpy.linspace(0, self.step_interval, self.step_number)
        self.electron_temperature = self.interpolate_profile(self.profiles.electron_temperature_2d)
        self.ion_temperature = self.interpolate_profile(self.profiles.ion_temperature_2d)
        self.density = self.interpolate_profile(self.profiles.density_2d)
        self.time_interval = [0, 0.005]
        self.time_slices = numpy.array([self.profiles.t_axis[time] for time in range(len(self.profiles.t_axis))
                                       if self.time_interval[0] < self.profiles.t_axis[time] < self.time_interval[1]])

    def interpolate_profile(self, profile):
        new_profile = interpn((self.profiles.z_axis, self.profiles.r_axis), profile, (self.vertical, self.steps))
        return new_profile

class Inputs1:
    initial_condition = numpy.array([1, 0, 0, 0, 0, 0, 0, 0, 0])
    step_interval = 0.1
    step_number = 100
    beam_energy = 60
    beam_species = 'Li'
    electron_temperature = numpy.linspace(1,6,step_number)*1e3
    proton_temperature = numpy.linspace(1,6,step_number)*1e3
    density = numpy.linspace(1,6,step_number)*1e20

    number_of_levels = len(initial_condition)
    steps = numpy.linspace(0, step_interval, step_number)


class Profiles:
    def __init__(self, filename, time_index):
        self.time_index = time_index
        self.filename = filename
        self.r_axis = numpy.array([0])
        self.z_axis = numpy.array([0])
        self.t_axis = numpy.array([0])
        self.density_2d = numpy.array([0])
        self.electron_temperature_2d = numpy.array([0])
        self.ion_temperature_2d = numpy.array([0])
        self.get_grid()
        self.get_density_2d(time_index=self.time_index)
        self.get_electron_temperature_2d(time_index=self.time_index)
        self.get_ion_temperature_2d(time_index=self.time_index)

    def get_grid(self):
        self.r_axis = get_data_from_hdf5.get_data_from_hdf5(self.filename, 'grid/xAxis')
        self.z_axis = get_data_from_hdf5.get_data_from_hdf5(self.filename, 'grid/yAxis')
        self.t_axis = get_data_from_hdf5.get_data_from_hdf5(self.filename, 'grid/tAxis')

    def get_density_2d(self,time_index):
        time_density_2d = get_data_from_hdf5.get_data_from_hdf5(self.filename, 'fields/density')
        self.density_2d = time_density_2d[time_index, :, :] # [t, z, r]
        del time_density_2d

    def get_electron_temperature_2d(self, time_index):
        time_temperature_2d = get_data_from_hdf5.get_data_from_hdf5(self.filename, 'fields/electronTemperature')
        self.electron_temperature_2d = time_temperature_2d[time_index, :, :]
        del time_temperature_2d

    def get_ion_temperature_2d(self, time_index):
        time_temperature_2d = get_data_from_hdf5.get_data_from_hdf5(self.filename, 'fields/ionTemperature')
        self.ion_temperature_2d = time_temperature_2d[time_index, :, :]
        del time_temperature_2d

    def calculate_new_r_axis(self):
        r = self.r_axis
        shifted_r = abs(r - max(r))
        self.r_axis = shifted_r[::-1]