Merge pull request #38 from gergopokol/development

Development merge into Master, ready for release 1.0.0.

.gitignore

@@ -1,4 +1,6 @@
+### PyCharm specific ###
 .idea/
+venv/
 
 ### Linux ###
 *~

.travis.yml

@@ -0,0 +1,7 @@
+language: python
+python:
+  - "3.6.1"
+install:
+  - pip install -r requirements.txt
+script:
+  - python -m unittest -v crm_solver.odetest.OdeTest

beamlet/benchmark_scenarios_h5.py

@@ -0,0 +1,21 @@
+import numpy
+import pandas
+
+
+class Scenarios:
+    def __init__(self):
+        self.grid = numpy.linspace(0,1,101)
+
+    def make_homo(self, component, q, z, a, density, temperature):
+        elements = pandas.DataFrame([[-1, q], [0, z], [0, a]], index=['q', 'Z', 'A'], columns=['electron', 'ion1']).transpose()
+
+        new_profiles = pandas.DataFrame([self.grid, [density]*len(self.grid), [temperature]*len(self.grid),
+                                         [density]*len(self.grid), [temperature]*len(self.grid)],
+                                        pandas.MultiIndex.from_arrays([['beamlet_grid', 'electron', 'electron', 'ion1', 'ion1'],
+                                                                       ['', 'density', 'temperature', 'density', 'temperature'],
+                                                                       ['m', 'm-3', 'eV', 'm-3', 'eV']],names=['type','property', 'unit'])).transpose()
+
+        elements.to_hdf('data/output/'+component+'_'+str(density)+'m-3_'+str(temperature)+'eV.h5', key='components')
+        new_profiles.to_hdf('data/output/'+component+'_'+str(density)+'m-3_'+str(temperature)+'eV.h5', key='profiles')
+
+

crm_solver/beamlet.py

@@ -8,12 +8,13 @@
 
 
 class Beamlet:
-    def __init__(self, param=None, profiles=None, data_path="beamlet/test0001.xml"):
+    def __init__(self, param=None, profiles=None, components=None, data_path="beamlet/testimp0001.xml"):
         self.param = param
         if not isinstance(self.param, etree._ElementTree):
             self.read_beamlet_param(data_path)
         self.profiles = profiles
-        if not isinstance(self.profiles, pandas.DataFrame):
+        self.components = components
+        if not (isinstance(self.components, pandas.DataFrame) and isinstance(self.profiles, pandas.DataFrame)):
             self.read_beamlet_profiles()
         if not isinstance(self.param.getroot().find('body').find('beamlet_mass'), etree._Element):
             self.get_mass()
@@ -28,10 +29,13 @@ def read_beamlet_param(self, data_path):
         print('Beamlet.param read from file: ' + data_path)
 
     def read_beamlet_profiles(self):
-        hdf5_path = self.param.getroot().find('body').find('beamlet_profiles').text
+        hdf5_path = self.param.getroot().find('body').find('beamlet_source').text
+        self.components = utility.getdata.GetData(data_path_name=hdf5_path, data_key=['components']).data
+        assert isinstance(self.components, pandas.DataFrame)
+        print('Beamlet.imp_components read from file: ' + hdf5_path)
         self.profiles = utility.getdata.GetData(data_path_name=hdf5_path, data_key=['profiles']).data
         assert isinstance(self.profiles, pandas.DataFrame)
-        print('Beamlet.profiles read from file: ' + hdf5_path)
+        print('Beamlet.imp_profiles read from file: ' + hdf5_path)
 
     def get_mass(self):
         data_path_name = 'atomic_data/' + self.param.getroot().find('body').find('beamlet_species').text + \
@@ -60,17 +64,17 @@ def get_velocity(self):
         return
 
     def initialize_ode(self):
-        self.coefficient_matrix = CoefficientMatrix(self.param, self.profiles)
+        self.coefficient_matrix = CoefficientMatrix(self.param, self.profiles, self.components)
         self.initial_condition = [1] + [0] * (self.coefficient_matrix.number_of_levels - 1)
 
     def solve_numerically(self):
         if self.coefficient_matrix is None or self.initial_condition is None:
             self.initialize_ode()
-        ode = Ode(coefficient_matrix=self.coefficient_matrix.matrix, initial_condition=self.initial_condition,
-                  steps=self.profiles['beamlet_grid'])
-        numerical = ode.calculate_solution()
+
+        ode = Ode(coeff_matrix=self.coefficient_matrix.matrix, init_condition=self.initial_condition)
+        numerical = ode.calculate_numerical_solution(self.profiles['beamlet grid']['distance']['m'])
+
         for level in range(self.coefficient_matrix.number_of_levels):
             label = 'level ' + str(level)
             self.profiles[label] = numerical[:, level]
         return
-

crm_solver/coefficientmatrix.py

@@ -3,18 +3,19 @@
 
 
 class CoefficientMatrix:
-    def __init__(self, beamlet_param, beamlet_profiles):
+    def __init__(self, beamlet_param, beamlet_profiles, plasma_components):
         self.beamlet_profiles = beamlet_profiles
-        self.rates = Rates(beamlet_param, beamlet_profiles)
+        self.rates = Rates(beamlet_param, beamlet_profiles, plasma_components)
         self.number_of_steps = self.rates.number_of_steps
         self.number_of_levels = self.rates.number_of_levels
         # Initialize matrices
         self.matrix = numpy.zeros(
             (self.number_of_levels, self.number_of_levels, self.number_of_steps))
         self.electron_terms = numpy.zeros(
             (self.number_of_levels, self.number_of_levels, self.number_of_steps))
-        self.ion_terms = numpy.zeros(
+        ion_terms = numpy.zeros(
             (self.number_of_levels, self.number_of_levels, self.number_of_steps))
+        self.ion_terms = numpy.concatenate([[ion_terms] * self.rates.number_of_ions])
         self.photon_terms = numpy.zeros(
             (self.number_of_levels, self.number_of_levels, self.number_of_steps))
         # Fill matrices
@@ -30,25 +31,28 @@ def assemble_terms(self):
                             - sum(self.rates.electron_neutral_collisions[from_level, :to_level, step]) \
                             - sum(self.rates.electron_neutral_collisions[from_level, (to_level+1):self.number_of_levels,
                                   step]) \
-                            - self.rates.electron_loss_collisions[0, from_level, step]
-                        self.ion_terms[from_level, to_level, step] = \
-                            - sum(self.rates.proton_neutral_collisions[from_level, :to_level, step]) \
-                            - sum(self.rates.proton_neutral_collisions[from_level, (to_level+1):self.number_of_levels,
-                                  step]) \
-                            - self.rates.electron_loss_collisions[1, from_level, step]
+                            - self.rates.electron_loss_collisions[from_level, step]
+                        for ion in range(self.rates.number_of_ions):
+                            self.ion_terms[ion][from_level, to_level, step] = \
+                                - sum(self.rates.ion_neutral_collisions[ion][from_level, :to_level, step]) \
+                                - sum(self.rates.ion_neutral_collisions[ion][from_level,
+                                      (to_level+1):self.number_of_levels, step]) \
+                                - self.rates.electron_loss_ion_collisions[ion][from_level, step]
                         self.photon_terms[from_level, to_level, step] = \
                             - sum(self.rates.einstein_coeffs[:, from_level]) / self.rates.velocity
                     else:
                         self.electron_terms[from_level, to_level, step] = \
                             self.rates.electron_neutral_collisions[from_level, to_level, step]
-                        self.ion_terms[from_level, to_level, step] = \
-                            self.rates.proton_neutral_collisions[from_level, to_level, step]
+                        for ion in range(self.rates.number_of_ions):
+                            self.ion_terms[ion][from_level, to_level, step] =\
+                                self.rates.ion_neutral_collisions[ion][from_level, to_level, step]
                         self.photon_terms[from_level, to_level, step] = \
                             self.rates.einstein_coeffs[to_level, from_level] / self.rates.velocity
 
     def apply_density(self):
         for step in range(self.number_of_steps):
-            self.matrix[:, :, step] = self.beamlet_profiles['beamlet_density'][step] * self.electron_terms[:, :, step] \
-                                      + self.beamlet_profiles['beamlet_density'][step] * self.ion_terms[:, :, step] \
-                                      + self.photon_terms[:, :, step]
-
+            self.matrix[:, :, step] = self.beamlet_profiles['electron']['density']['m-3'][step] * self.electron_terms[:, :, step]
+            for ion in range(self.rates.number_of_ions):
+                self.matrix[:, :, step] = self.matrix[:, :, step] + self.beamlet_profiles['ion' + str(ion+1)]['density']['m-3'][step] \
+                                                                    * self.ion_terms[ion][:, :, step]
+            self.matrix[:, :, step] = self.matrix[:, :, step] + self.photon_terms[:, :, step]

crm_solver/ode.py

@@ -4,39 +4,67 @@
 
 
 class Ode:
-    def __init__(self,
-                 coefficient_matrix=numpy.array([0.]),
-                 initial_condition=numpy.array([0.]),
-                 steps=numpy.array([0.])):
-        self.coefficient_matrix = coefficient_matrix
-        self.initial_condition = initial_condition
-        self.steps = steps
-
-    def calculate_solution(self):
-        solution = odeint(func=self.set_up_equation, y0=self.initial_condition, t=self.steps,
-                          args=(self.coefficient_matrix, self.steps))
-        return solution
-
-    def analytical_solution(self):
-        eigenvalues, eigenvectors = numpy.linalg.eig(self.coefficient_matrix)
-        if self.initial_condition.size == 1:
-            analytical_solution = numpy.zeros(self.steps.size)
-            for step in range(self.steps.size):
-                analytical_solution[step] = 1 / eigenvectors * self.initial_condition * eigenvectors \
-                                         * numpy.exp(eigenvalues * self.steps[step])
+    def __init__(self, coeff_matrix, init_condition):
+        self.coeff_matrix = coeff_matrix
+        self.init_condition = init_condition
+        self.__validate_parameters()
+
+    def __validate_parameters(self):
+        self.__validate(self.coeff_matrix)
+        self.__validate(self.init_condition)
+        self.__validate_matrix(self.coeff_matrix)
+
+    @staticmethod
+    def __validate(parameter):
+        if parameter is None:
+            raise ValueError('Try to define Ode class with null matrix.')
+
+    def __validate_matrix(self, parameter):
+        if self.__is_not_square(parameter):
+            raise ValueError("The matrix must be squared")
+
+    @staticmethod
+    def __is_not_square(matrix):
+        matrix_temp1 = matrix.shape[0]
+        if matrix.ndim > 1:
+            matrix_temp2 = matrix.shape[1]
         else:
-            analytical_solution = numpy.zeros((self.steps.size, self.initial_condition.size))
-            for step in range(self.steps.size):
-                analytical_solution[step, :] = numpy.dot(numpy.dot(numpy.linalg.inv(eigenvectors),
-                                                                   self.initial_condition), eigenvectors)\
-                                               * numpy.exp(eigenvalues * self.steps[step])
+            matrix_temp2 = 0
+        return matrix_temp1 != matrix_temp2
+
+    def calculate_numerical_solution(self, steps):
+        return odeint(func=self.set_derivative_vector, y0=self.init_condition, t=steps, args=(self.coeff_matrix, steps))
+
+    def calculate_analytical_solution(self, steps):
+        eigenvalues, eigenvectors = numpy.linalg.eig(self.coeff_matrix)
+        if self.init_condition.size == 1:
+            return self.__calculate_analytical_solution_1d(steps, eigenvalues)
+        else:
+            return self.__calculate_analytical_solution_else(steps, eigenvalues, eigenvectors)
+
+    def __calculate_analytical_solution_else(self, steps, eigenvalues, eigenvectors):
+        analytical_solution = numpy.zeros((steps.size, self.init_condition.size))
+        temporal_constant = numpy.dot(numpy.dot(numpy.linalg.inv(eigenvectors), self.init_condition), eigenvectors)
+        for step in range(steps.size):
+            analytical_solution[step, :] = self.calculate_exp_solution(temporal_constant, eigenvalues, steps[step])
+        return analytical_solution
+
+    def __calculate_analytical_solution_1d(self, steps, eigenvalues):
+        analytical_solution = numpy.zeros(steps.size)
+        for step in range(steps.size):
+            analytical_solution[step] = self.calculate_exp_solution(self.init_condition, eigenvalues, steps[step])
         return analytical_solution
 
     @staticmethod
-    def set_up_equation(variable_vector, actual_position, coefficient_matrix, steps):
-        if coefficient_matrix.ndim == 3:
-            interp_coefficient_matrix = interp1d(steps, coefficient_matrix, axis=2, fill_value='extrapolate')
-            derivative_vector = numpy.dot(variable_vector, interp_coefficient_matrix(actual_position))
-        elif coefficient_matrix.ndim == 2:
-            derivative_vector = numpy.dot(variable_vector, coefficient_matrix)
-        return derivative_vector
+    def set_derivative_vector(variable_vector, actual_position, coeff_matrix, steps):
+        if coeff_matrix.ndim == 3:
+            interp_coeff_matrix = interp1d(steps, coeff_matrix, axis=2, fill_value='extrapolate')
+            return numpy.dot(variable_vector, interp_coeff_matrix(actual_position))
+        elif coeff_matrix.ndim == 2:
+            return numpy.dot(variable_vector, coeff_matrix)
+        else:
+            raise ValueError('Rate Coefficient Matrix of dimensions: ' + str(coeff_matrix.ndim) + ' is not supported')
+
+    @staticmethod
+    def calculate_exp_solution(init, coefficient, variable):
+        return init * numpy.exp(coefficient * variable)