Perturbative Raman Solver

Raman effect evaluation based on a perturbative solution for faster computation

Change-Id: Ie6d4ea4d9f95d8755dc8dfd004c954d4c2c5f759

docs/json.rst

@@ -492,6 +492,19 @@ See ``delta_power_range_db`` for more explaination.
 |                                             |           | mandatory when Raman amplification is       |
 |                                             |           | included in the simulation                  |
 +---------------------------------------------+-----------+---------------------------------------------+
+| ``raman_params.method``                     | (string)  | Model used for Raman evaluation. Valid      |
+|                                             |           | choices are ``perturbative`` (see           |
+|                                             |           | `arXiv:2304.11756                           |
+|                                             |           | <https://arxiv.org/abs/2304.11756>`_) and   |
+|                                             |           | ``numerical``, the GNPy legacy first order  |
+|                                             |           | derivative numerical solution.              |
++---------------------------------------------+-----------+---------------------------------------------+
+|``raman_params.order``                       |           | Order of the perturbative expansion.        |
+|                                             |           | For C- and C+L-band transmission scenarios  |
+|                                             |           | the second order provides high accuracy     |
+|                                             |           | considering common values of fiber input    |
+|                                             |           | power. (Default is 2)                       |
++---------------------------------------------+-----------+---------------------------------------------+
 | ``raman_params.result_spatial_resolution``  | (number)  | Spatial resolution of the output            |
 |                                             |           | Raman profile along the entire fiber span.  |
 |                                             |           | This affects the accuracy and the           |
@@ -503,11 +516,18 @@ See ``delta_power_range_db`` for more explaination.
 |                                             |           | channel around 0 dBm, a suggested value of  |
 |                                             |           | spatial resolution is 10e3 m                |
 +---------------------------------------------+-----------+---------------------------------------------+
-| ``raman_params.solver_spatial_resolution``  | (number)  | Spatial step for the iterative solution     |
-|                                             |           | of the first order differential equation    |
-|                                             |           | used to calculate the Raman profile         |
-|                                             |           | along the entire fiber span.                |
-|                                             |           | This affects the accuracy and the           |
+| ``raman_params.solver_spatial_resolution``  | (number)  | When using the ``perturbative`` method,     |
+|                                             |           | the step for the spatial integration does   |
+|                                             |           | not affect the first order. Therefore, a    |
+|                                             |           | large step can be used when no              |
+|                                             |           | counter-propagating Raman amplification is  |
+|                                             |           | present; a suggested value is 10e3 m.       |
+|                                             |           | In presence of counter-propagating Raman    |
+|                                             |           | amplification or when using the             |
+|                                             |           | ``numerical`` method the following remains  |
+|                                             |           | valid.                                      |
+|                                             |           | The spatial step for the iterative solution |
+|                                             |           | affects the accuracy and the                |
 |                                             |           | computational time of the evaluated         |
 |                                             |           | Raman profile:                              |
 |                                             |           | smaller the spatial resolution higher both  |

gnpy/core/parameters.py

@@ -36,25 +36,32 @@ def __init__(self, power, frequency, propagation_direction):
 
 
 class RamanParams(Parameters):
-    def __init__(self, flag=False, result_spatial_resolution=10e3, solver_spatial_resolution=50):
+    def __init__(self, flag=False, method='perturbative', order=2, result_spatial_resolution=10e3,
+                 solver_spatial_resolution=10e3):
         """Simulation parameters used within the Raman Solver
 
         :params flag: boolean for enabling/disable the evaluation of the Raman power profile in frequency and position
+        :params method: Raman solver method
+        :params order: solution order for perturbative method
         :params result_spatial_resolution: spatial resolution of the evaluated Raman power profile
         :params solver_spatial_resolution: spatial step for the iterative solution of the first order ode
         """
         self.flag = flag
+        self.method = method
+        self.order = order
         self.result_spatial_resolution = result_spatial_resolution  # [m]
         self.solver_spatial_resolution = solver_spatial_resolution  # [m]
 
     def to_json(self):
         return {"flag": self.flag,
+                "method": self.method,
+                "order": self.order,
                 "result_spatial_resolution": self.result_spatial_resolution,
                 "solver_spatial_resolution": self.solver_spatial_resolution}
 
 
 class NLIParams(Parameters):
-    def __init__(self, method='gn_model_analytic', dispersion_tolerance=1, phase_shift_tolerance=0.1,
+    def __init__(self, method='gn_model_analytic', dispersion_tolerance=4, phase_shift_tolerance=0.1,
                  computed_channels=None, computed_number_of_channels=None):
         """Simulation parameters used within the Nli Solver
 

gnpy/core/science_utils.py

@@ -15,11 +15,12 @@
     polyfit, log, reshape, swapaxes, full, nan
 from logging import getLogger
 from scipy.constants import k, h
+from scipy.integrate import cumtrapz
 from scipy.interpolate import interp1d
-from math import isclose
+from math import isclose, factorial
 
 from gnpy.core.utils import db2lin, lin2db
-from gnpy.core.exceptions import EquipmentConfigError
+from gnpy.core.exceptions import EquipmentConfigError, ParametersError
 from gnpy.core.parameters import SimParams
 from gnpy.core.info import SpectralInformation
 
@@ -136,15 +137,17 @@ def calculate_stimulated_raman_scattering(spectral_info: SpectralInformation, fi
                     co_cr = fiber.cr(co_frequency)
                     co_alpha = fiber.alpha(co_frequency)
                     co_power_profile = \
-                        RamanSolver.first_order_derivative_solution(co_power, co_alpha, co_cr, z, lumped_losses)
+                        RamanSolver.calculate_unidirectional_stimulated_raman_scattering(co_power, co_alpha, co_cr, z,
+                                                                                         lumped_losses)
                 # Counter-propagating profile initialization
                 cnt_power_profile = zeros([cnt_frequency.size, z.size])
                 if cnt_frequency.size:
                     cnt_cr = fiber.cr(cnt_frequency)
                     cnt_alpha = fiber.alpha(cnt_frequency)
-                    cnt_power_profile = \
-                        flip(RamanSolver.first_order_derivative_solution(cnt_power, cnt_alpha, cnt_cr,
-                                                                         z[-1] - flip(z), flip(lumped_losses)), axis=1)
+                    cnt_power_profile = flip(
+                        RamanSolver.calculate_unidirectional_stimulated_raman_scattering(cnt_power, cnt_alpha, cnt_cr,
+                                                                                         z[-1] - flip(z),
+                                                                                         flip(lumped_losses)), axis=1)
                 # Co-propagating and Counter-propagating Profile Computation
                 if co_frequency.size and cnt_frequency.size:
                     co_power_profile, cnt_power_profile = \
@@ -163,8 +166,9 @@ def calculate_stimulated_raman_scattering(spectral_info: SpectralInformation, fi
                 alpha = fiber.alpha(spectral_info.frequency)
                 cr = fiber.cr(spectral_info.frequency)
                 # Power profile
-                power_profile = \
-                    RamanSolver.first_order_derivative_solution(spectral_info.signal, alpha, cr, z, lumped_losses)
+                power_profile = (
+                    RamanSolver.calculate_unidirectional_stimulated_raman_scattering(spectral_info.signal, alpha, cr, z,
+                                                                                     lumped_losses))
                 # Loss profile
                 loss_profile = power_profile / outer(spectral_info.signal, ones(z.size))
                 frequency = spectral_info.frequency
@@ -200,8 +204,8 @@ def calculate_spontaneous_raman_scattering(spectral_info: SpectralInformation, s
         return ase
 
     @staticmethod
-    def first_order_derivative_solution(power_in, alpha, cr, z, lumped_losses):
-        """Solves the Raman first order derivative equation
+    def calculate_unidirectional_stimulated_raman_scattering(power_in, alpha, cr, z, lumped_losses):
+        """Solves the Raman equation
 
         :param power_in: launch power array
         :param alpha: loss coefficient array
@@ -210,18 +214,58 @@ def first_order_derivative_solution(power_in, alpha, cr, z, lumped_losses):
         :param lumped_losses: concentrated losses array along the fiber span
         :return: power profile matrix
         """
-        dz = z[1:] - z[:-1]
-        power = outer(power_in, ones(z.size))
-        for i in range(1, z.size):
-            power[:, i] = \
-                power[:, i - 1] * (1 + (- alpha + sum(cr * power[:, i - 1], 1)) * dz[i - 1]) * lumped_losses[i - 1]
+        if sim_params.raman_params.method == 'perturbative':
+            if sim_params.raman_params.order > 4:
+                raise ValueError(f'Order {sim_params.raman_params.order} not implemented in Raman Solver.')
+            z_lumped_losses = append(z[lumped_losses != 1], z[-1])
+            llumped_losses = append(1, lumped_losses[lumped_losses != 1])
+            power = outer(power_in, ones(z.size))
+            last_position = 0
+            z_indices = arange(0, z.size)
+
+            for z_lumped_loss, lumped_loss in zip(z_lumped_losses, llumped_losses):
+                if last_position < z[-1]:
+                    interval = z_indices[(z >= last_position) * (z <= z_lumped_loss) == 1]
+                    z_interval = z[interval] - last_position
+                    last_position = z[interval][-1]
+                    p0 = power_in * lumped_loss
+                    power_interval = outer(p0, ones(z_interval.size))
+                    alphaz = outer(alpha, z_interval)
+                    expz = exp(- alphaz)
+                    eff_length = 1 / outer(alpha, ones(z_interval.size)) * (1 - expz)
+                    crpz = transpose(ones([z_interval.size, cr.shape[0], cr.shape[1]]) * cr * p0, (1, 2, 0))
+                    exponent = - alphaz
+                    if sim_params.raman_params.order >= 1:
+                        gamma1 = sum(crpz * eff_length, 1)
+                        exponent += gamma1
+                    if sim_params.raman_params.order >= 2:
+                        gamma2 = sum(crpz * cumtrapz(expz * gamma1, z_interval, axis=1, initial=0), 1)
+                        exponent += gamma2
+                    if sim_params.raman_params.order >= 3:
+                        gamma3 = sum(crpz * cumtrapz(expz * (gamma2 + 1/2 * gamma1**2), z_interval,
+                                                     axis=1, initial=0), 1)
+                        exponent += gamma3
+                    if sim_params.raman_params.order >= 4:
+                        gamma4 = sum(crpz * cumtrapz(expz * (gamma3 + gamma1 * gamma2 + 1/factorial(3) * gamma1**3),
+                                                     z_interval, axis=1, initial=0), 1)
+                        exponent += gamma4
+                    power_interval *= exp(exponent)
+                    power[:, interval[1:]] = power_interval[:, 1:]
+                    power_in = power_interval[:, -1]
+        elif sim_params.raman_params.method == 'numerical':
+            dz = z[1:] - z[:-1]
+            power = outer(power_in, ones(z.size))
+            for i in range(1, z.size):
+                power[:, i] = (power[:, i - 1] * (1 + (- alpha + sum(cr * power[:, i - 1], 1)) * dz[i - 1]) *
+                               lumped_losses[i - 1])
+        else:
+            raise ValueError(f'Method {sim_params.raman_params.method} not implemented in Raman Solver.')
         return power
 
     @staticmethod
     def iterative_algorithm(co_initial_guess_power, cnt_initial_guess_power, co_frequency, cnt_frequency, z, fiber,
                             lumped_losses):
-        """Solves the Raman first order derivative equation in case of both co- and counter-propagating
-        frequencies
+        """Solves the Raman equation in case of both co- and counter-propagating frequencies
 
         :param co_initial_guess_power: co-propagationg Raman first order derivative equation solution
         :param cnt_initial_guess_power: counter-propagationg Raman first order derivative equation solution

tests/data/sim_params.json

@@ -2,7 +2,7 @@
   "raman_params": {
     "flag": true,
     "result_spatial_resolution": 10e3,
-    "solver_spatial_resolution": 5
+    "solver_spatial_resolution": 10e3
   },
   "nli_params": {
     "method": "ggn_spectrally_separated",

tests/test_network_functions.py

@@ -553,13 +553,13 @@ def test_get_node_restrictions(cls, defaultparams, variety_list, booster_list, b
     ('no_design', 'Multiband_amplifier', 'LBAND', 10.0, 0.0, 'std_low_gain_multiband_bis', False),
     ('type_variety', 'Multiband_amplifier', 'LBAND', 10.0, 0.0, 'std_medium_gain_multiband', False),
     ('design', 'Multiband_amplifier', 'LBAND', 9.344985, 0.0, 'std_medium_gain_multiband', True),
-    ('no_design', 'Multiband_amplifier', 'LBAND', 9.344985, -0.94256, 'std_low_gain_multiband_bis', True),
-    ('no_design', 'Multiband_amplifier', 'CBAND', 10.980212, -1.60348, 'std_low_gain_multiband_bis', True),
+    ('no_design', 'Multiband_amplifier', 'LBAND', 9.344985, -0.938676, 'std_low_gain_multiband_bis', True),
+    ('no_design', 'Multiband_amplifier', 'CBAND', 10.977065, -1.600193, 'std_low_gain_multiband_bis', True),
     ('no_design', 'Fused', 'LBAND', 21.0, 0.0, 'std_medium_gain_multiband', False),
-    ('no_design', 'Fused', 'LBAND', 20.344985, -0.82184, 'std_medium_gain_multiband', True),
-    ('no_design', 'Fused', 'CBAND', 21.773072, -1.40300, 'std_medium_gain_multiband', True),
-    ('design', 'Fused', 'CBAND', 21.214048, 0.0, 'std_medium_gain_multiband', True),
-    ('design', 'Multiband_amplifier', 'CBAND', 11.044233, 0.0, 'std_medium_gain_multiband', True)])
+    ('no_design', 'Fused', 'LBAND', 20.344985, -0.819176, 'std_medium_gain_multiband', True),
+    ('no_design', 'Fused', 'CBAND', 21.770319, -1.40032, 'std_medium_gain_multiband', True),
+    ('design', 'Fused', 'CBAND', 21.21108, 0.0, 'std_medium_gain_multiband', True),
+    ('design', 'Multiband_amplifier', 'CBAND', 11.041037, 0.0, 'std_medium_gain_multiband', True)])
 def test_multiband(case, site_type, band, expected_gain, expected_tilt, expected_variety, sim_params):
     """Check:
     - if amplifiers are defined in multiband they are used for design,
@@ -612,8 +612,10 @@ def test_tilt_fused():
         estimate_srs_power_deviation(network, node, equipment, design_bands, input_powers)
     # restore simParams
     SimParams.set_params(save_sim_params)
-    assert fused_tilt_db == tilt_db
-    assert fused_tilt_target == tilt_target
+    for key in tilt_db:
+        assert_allclose(tilt_db[key], fused_tilt_db[key], rtol=1e-3)
+    for key in tilt_target:
+        assert_allclose(tilt_target[key], fused_tilt_target[key], rtol=1e-3)
 
 
 def network_wo_booster(site_type, bands):

tests/test_science_utils.py

@@ -63,7 +63,7 @@ def test_fiber():
     assert_allclose(p_nli, expected_results['nli'], rtol=1e-3)
 
     # propagation with Raman
-    SimParams.set_params({'raman_params': {'flag': True, 'solver_spatial_resolution': 1}})
+    SimParams.set_params({'raman_params': {'flag': True}})
     spectral_info_out = fiber(spectral_info_input)
 
     p_signal = spectral_info_out.signal
@@ -82,6 +82,7 @@ def test_raman_fiber():
                                                             baud_rate=32e9, spacing=50e9, tx_osnr=40.0,
                                                             tx_power=1e-3)
     SimParams.set_params(load_json(TEST_DIR / 'data' / 'sim_params.json'))
+    SimParams().raman_params.solver_spatial_resolution = 5
     fiber = RamanFiber(**load_json(TEST_DIR / 'data' / 'test_science_utils_fiber_config.json'))
     fiber.ref_pch_in_dbm = 0.0
     # propagation
@@ -120,18 +121,19 @@ def test_fiber_lumped_losses_srs(set_sim_params):
                                                             baud_rate=32e9, spacing=50e9, tx_osnr=40.0,
                                                             tx_power=1e-3)
 
-    SimParams.set_params(load_json(TEST_DIR / 'data' / 'sim_params.json'))
     fiber = Fiber(**load_json(TEST_DIR / 'data' / 'test_lumped_losses_raman_fiber_config.json'))
     raman_fiber = RamanFiber(**load_json(TEST_DIR / 'data' / 'test_lumped_losses_raman_fiber_config.json'))
 
     # propagation
     # without Raman pumps
+    SimParams.set_params({'raman_params': {'flag': True, 'order': 2}})
     stimulated_raman_scattering = RamanSolver.calculate_stimulated_raman_scattering(spectral_info_input, fiber)
     power_profile = stimulated_raman_scattering.power_profile
     expected_power_profile = read_csv(TEST_DIR / 'data' / 'test_lumped_losses_fiber_no_pumps.csv').values
     assert_allclose(power_profile, expected_power_profile, rtol=1e-3)
 
     # with Raman pumps
+    SimParams.set_params({'raman_params': {'flag': True, 'order': 1, 'solver_spatial_resolution': 5}})
     stimulated_raman_scattering = RamanSolver.calculate_stimulated_raman_scattering(spectral_info_input, raman_fiber)
     power_profile = stimulated_raman_scattering.power_profile
     expected_power_profile = read_csv(TEST_DIR / 'data' / 'test_lumped_losses_raman_fiber.csv').values