Squashed commit of the following:

commit 377d0f5 Merge: 3448172 e4bc7ec Author: retkiewi <58555777+retkiewi@users.noreply.github.com> Date: Sun Jun 9 15:50:22 2024 +0200 Merge pull request #55 from jbrage/feature/43-reorganize-hadrons Feature/43 reorganize electrons commit e4bc7ec Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat Jun 8 19:57:37 2024 +0200 improve runners, beam simulation and use more numpy commit 532ac8e Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 25 20:45:01 2024 +0200 fix nit commit 13ed258 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 25 20:39:26 2024 +0200 move calculate method to generic solver commit 08bc588 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 25 20:31:17 2024 +0200 only return f in run simulation for continous commit f2157e3 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 25 20:22:43 2024 +0200 unify calculate method for pulsed and continous beams commit 9a1247c Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 25 20:10:25 2024 +0200 return f as an array of results after each time step commit 2323a81 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 25 19:27:42 2024 +0200 fix getting solver name in runner commit 5ed8619 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 25 19:14:38 2024 +0200 update simulation runners commit a2c3e3b Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 25 18:32:33 2024 +0200 use dataclasses for solvers commit 55aa3aa Author: retkiewi <krzret@student.agh.edu.pl> Date: Sun May 19 16:32:29 2024 +0200 add simulation runners, fix bugs commit 3240b58 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sun May 19 15:29:12 2024 +0200 optimize constant calculations commit 6c77bb7 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sun May 19 15:27:00 2024 +0200 remove parameter dict in generic solver commit 70a3ef5 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sun May 19 15:09:51 2024 +0200 fix generic_electron_solver formatting commit d29a577 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 18 19:28:05 2024 +0200 extract post beam electron density calculation to a separate method commit 1305350 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 18 19:18:20 2024 +0200 extract common init commit 841d56a Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 18 18:53:48 2024 +0200 small cleanup of pulsed_e_beam.py commit 68e9ff2 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 18 18:50:46 2024 +0200 convert electrons/continous_e_beam.py to python commit 6b4df88 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 18 18:42:14 2024 +0200 large cleanup in electrons/continous_e_beam.py commit 541d5cc Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 18 18:41:46 2024 +0200 small cleanup in initial_recombination.py commit a404c04 Author: retkiewi <krzret@student.agh.edu.pl> Date: Sat May 18 18:03:00 2024 +0200 convert electrons/continous_e_beam.py to python commit 3448172 Author: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon Apr 1 15:23:06 2024 +0000 Bump actions/setup-python from 5.0.0 to 5.1.0 Bumps [actions/setup-python](https://github.com/actions/setup-python) from 5.0.0 to 5.1.0. - [Release notes](https://github.com/actions/setup-python/releases) - [Commits](actions/setup-python@v5.0.0...v5.1.0) --- updated-dependencies: - dependency-name: actions/setup-python dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] <support@github.com> commit d6c0035 Author: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon Dec 11 15:54:38 2023 +0000 Bump actions/setup-python from 4.7.1 to 5.0.0 Bumps [actions/setup-python](https://github.com/actions/setup-python) from 4.7.1 to 5.0.0. - [Release notes](https://github.com/actions/setup-python/releases) - [Commits](actions/setup-python@v4.7.1...v5.0.0) --- updated-dependencies: - dependency-name: actions/setup-python dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] <support@github.com> commit a246282 Author: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon Oct 2 20:37:55 2023 +0000 Bump actions/checkout from 3 to 4 Bumps [actions/checkout](https://github.com/actions/checkout) from 3 to 4. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](actions/checkout@v3...v4) --- updated-dependencies: - dependency-name: actions/checkout dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] <support@github.com> commit 8ef4655 Author: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon Oct 2 15:30:14 2023 +0000 Bump actions/setup-python from 4.5.0 to 4.7.1 Bumps [actions/setup-python](https://github.com/actions/setup-python) from 4.5.0 to 4.7.1. - [Release notes](https://github.com/actions/setup-python/releases) - [Commits](actions/setup-python@v4.5.0...v4.7.1) --- updated-dependencies: - dependency-name: actions/setup-python dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] <support@github.com> commit 0dce892 Author: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon Jan 16 15:46:26 2023 +0000 Bump actions/setup-python from 4.4.0 to 4.5.0 Bumps [actions/setup-python](https://github.com/actions/setup-python) from 4.4.0 to 4.5.0. - [Release notes](https://github.com/actions/setup-python/releases) - [Commits](actions/setup-python@v4.4.0...v4.5.0) --- updated-dependencies: - dependency-name: actions/setup-python dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] <support@github.com> commit b5ca6d8 Author: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon Dec 26 15:06:08 2022 +0000 Bump actions/setup-python from 4.3.1 to 4.4.0 Bumps [actions/setup-python](https://github.com/actions/setup-python) from 4.3.1 to 4.4.0. - [Release notes](https://github.com/actions/setup-python/releases) - [Commits](actions/setup-python@v4.3.1...v4.4.0) --- updated-dependencies: - dependency-name: actions/setup-python dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] <support@github.com> commit 24ecd26 Author: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon Dec 12 15:08:26 2022 +0000 Bump actions/setup-python from 4.3.0 to 4.3.1 Bumps [actions/setup-python](https://github.com/actions/setup-python) from 4.3.0 to 4.3.1. - [Release notes](https://github.com/actions/setup-python/releases) - [Commits](actions/setup-python@v4.3.0...v4.3.1) --- updated-dependencies: - dependency-name: actions/setup-python dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] <support@github.com> commit d2f22ba Merge: ee90579 70d9bce Author: retkiewi <58555777+retkiewi@users.noreply.github.com> Date: Thu Dec 8 20:45:07 2022 +0100 Merge pull request #45 from jbrage/feature/43-reorganize-hadrons Reorganize hadrons package
jbrage · Jun 14, 2024 · 1c4e36e · 1c4e36e
1 parent 3656864
commit 1c4e36e
Show file tree

Hide file tree

Showing 8 changed files with 455 additions and 15 deletions.
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
@@ -21,10 +21,10 @@ jobs:
          python-version: ['3.8', '3.9', '3.10']
          platform: [ubuntu-latest]  # TODO add Windows as well
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - name: Set up Python
-        uses: actions/setup-python@v4.3.0
+        uses: actions/setup-python@v5.1.0
         with:
           python-version: ${{ matrix.python-version }}
 

diff --git a/electrons/cython/run_simulation.py b/electrons/cython/run_simulation.py
@@ -0,0 +1,69 @@
+from electrons.cython.pulsed_e_beam import pulsed_beam_PDEsolver
+from electrons.cython.continuous_e_beam import continuous_beam_PDEsolver
+
+import argparse
+
+SOLVER_MAP = {"continous": continuous_beam_PDEsolver, "pulsed": pulsed_beam_PDEsolver}
+
+
+def run_simulation(
+    solver_name="continous",
+    voltage_V=300,
+    electrode_gap=0.1,
+    electron_density_per_cm3=1e9,
+    verbose=True,
+):
+    parameters = {
+        "voltage_V": voltage_V,
+        "d_cm": electrode_gap,
+        "elec_per_cm3": electron_density_per_cm3,
+        "show_plot": False,
+        "print_parameters": False,
+    }
+
+    if solver_name not in SOLVER_MAP.keys():
+        print(f'Invalid solver type "{solver_name}", defaulting to Continous solver.')
+        solver_name = "continous"
+
+    if verbose:
+        print(f"Running the simulation using the {solver_name} solver.")
+        print(f"Voltage: {voltage_V} [V]")
+        print(f"Electrode gap: {electrode_gap} [cm]")
+        print(f"Electron density per cm3: {electron_density_per_cm3}")
+
+    solver = SOLVER_MAP[solver_name]
+
+    # return the collection efficiency
+    result = solver(parameters)
+
+    if solver_name == "continous":
+        f = result[1]["f"]
+    else:
+        f = result
+
+    return f
+
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument(
+        "solver_name",
+        type=str,
+        default="continous",
+        help="The type of the solver to use",
+    )
+
+    parser.add_argument(
+        "--verbose",
+        "-v",
+        type=bool,
+        default=True,
+    )
+
+    args = parser.parse_args()
+
+    result = run_simulation(**vars(args))
+
+    print("calculated f is", result)
diff --git a/electrons/python/continuous_e_beam.py b/electrons/python/continuous_e_beam.py
@@ -0,0 +1,38 @@
+from __future__ import division
+
+from generic_electron_solver import GenericElectronSolver
+
+
+class ContinousBeamPDEsolver(GenericElectronSolver):
+    def get_initialised_charge_carriers_after_beam(self):
+        electron_density_per_cm3_s = self.electron_density_per_cm3 * self.dt
+
+        initialised_charge_carriers = (
+            electron_density_per_cm3_s
+            * (self.no_xy - 2 * self.delta_border) ** 2
+            * self.no_z
+        )
+
+        return initialised_charge_carriers
+
+    def update_electron_density_arrays_after_beam(
+        self,
+        positive_array,
+        negative_array,
+    ):
+        electron_density_per_cm3_s = self.electron_density_per_cm3 * self.dt
+
+        positive_array[
+            self.delta_border : (self.no_xy - self.delta_border),
+            self.delta_border : (self.no_xy - self.delta_border),
+            self.no_z_electrode : (self.no_z + self.no_z_electrode),
+        ] += electron_density_per_cm3_s
+
+        negative_array[
+            self.delta_border : (self.no_xy - self.delta_border),
+            self.delta_border : (self.no_xy - self.delta_border),
+            self.no_z_electrode : (self.no_z + self.no_z_electrode),
+        ] += electron_density_per_cm3_s
+
+    def should_simulate_beam_for_time_step(self, time_step):
+        return True
diff --git a/electrons/python/generic_electron_solver.py b/electrons/python/generic_electron_solver.py
@@ -0,0 +1,238 @@
+from __future__ import division
+from abc import ABC, abstractmethod
+from dataclasses import dataclass
+from numpy.typing import NDArray
+from typing import Tuple
+import numpy as np
+
+
+def von_neumann_expression(
+    dt, ion_diff, grid_spacing_cm, ion_mobility, Efield_V_cm
+) -> None:
+    """
+    Finds a time step dt which fulfils the von Neumann criterion, i.e. ensures the numericl error does not increase but
+    decreases and eventually damps out
+    """
+    von_neumann_expression = False
+
+    while not von_neumann_expression:
+        dt /= 1.01
+        # as defined in the Deghan (2004) paper
+        # sx = ion_diff*dt/(self.grid_spacing_cm**2)
+        # sy = ion_diff*dt/(self.grid_spacing_cm**2)
+        sx = sy = 0.0  # uniform charge => no gradient driven diffusion in the xy plane
+
+        sz = ion_diff * dt / (grid_spacing_cm**2)
+        cx = cy = 0.0
+        cz = ion_mobility * Efield_V_cm * dt / grid_spacing_cm
+        # check von Neumann's criterion
+        criterion_1 = 2 * (sx + sy + sz) + cx**2 + cy**2 + cz**2 <= 1
+
+        criterion_2 = cx**2 * cy**2 * cz**2 <= 8 * sx * sy * sz
+
+        von_neumann_expression = criterion_1 and criterion_2
+
+    return dt, sx, sy, sz, cx, cy, cz
+
+
+@dataclass
+class GenericElectronSolver(ABC):
+    # Simulation parameters
+    electron_density_per_cm3: float  # fluence-rate [/cm^2/s]
+    voltage_V: float  # [V/cm] magnitude of the electric field
+    electrode_gap: float  # [cm] # electrode gap
+    grid_spacing_cm: float
+    ion_mobility: float = 1.73  # cm s^-1 V^-1, averaged for positive and negative ions
+    ion_diff: float = 3.7e-2  # cm^2/s, averaged for positive and negative ions
+    alpha: float = 1.60e-6  # cm^3/s, recombination constant
+    r_cm: float = 0.002  # radius of the sampled cylinder
+    buffer_radius: int = 5
+    no_z_electrode: int = (
+        5  # length of the electrode-buffer to ensure no ions drift through the array in one time step
+    )
+    dt: float = 1.0
+
+    # Derived properties
+    @property
+    def no_xy(self) -> int:
+        """Number of voxels om the xy-directions"""
+
+        no_xy = int(2 * self.r_cm / self.grid_spacing_cm)
+        no_xy += 2 * self.buffer_radius
+        return no_xy
+
+    @property
+    def no_z(self) -> int:
+        """Number of voxels om the z-direction"""
+
+        return int(self.electrode_gap / self.grid_spacing_cm)
+
+    @property
+    def no_z_with_buffer(self) -> int:
+        return 2 * self.no_z_electrode + self.no_z
+
+    @property
+    def Efield_V_cm(self) -> float:
+        return self.voltage_V / self.electrode_gap
+
+    @property
+    def separation_time_steps(self) -> int:
+        """
+        Number of step required to drag the two charge carrier distributions apart
+        along with the number of tracks to be uniformly distributed over the domain
+        """
+        return int(
+            self.electrode_gap / (2.0 * self.ion_mobility * self.Efield_V_cm * self.dt)
+        )
+
+    @property
+    def computation_time_steps(self) -> int:
+        return self.separation_time_steps * 3
+
+    @property
+    def delta_border(self) -> int:
+        return 2
+
+    def __post_init__(
+        self,
+    ):
+        # depending on the cluster/computer, the upper limit may be changed
+        if (self.no_xy * self.no_xy * self.no_z) > 1e8:
+            raise ValueError(
+                "Too many elements in the array: %i"
+                % (self.no_xy * self.no_xy * self.no_z_with_buffer)
+            )
+
+        self.dt, self.sx, self.sy, self.sz, self.cx, self.cy, self.cz = (
+            von_neumann_expression(
+                self.dt,
+                self.ion_diff,
+                self.grid_spacing_cm,
+                self.ion_mobility,
+                self.Efield_V_cm,
+            )
+        )
+
+    @abstractmethod
+    def should_simulate_beam_for_time_step(self, time_step: int) -> bool:
+        pass
+
+    @abstractmethod
+    def update_electron_density_arrays_after_beam(
+        self, positive_array: NDArray, negative_array: NDArray
+    ) -> None:
+        pass
+
+    @abstractmethod
+    def get_initialised_charge_carriers_after_beam() -> float:
+        pass
+
+    def simulate_beam(self, positive_array: NDArray, negative_array: NDArray):
+        self.update_electron_density_arrays_after_beam(positive_array, negative_array)
+
+        return self.get_initialised_charge_carriers_after_beam()
+
+    def calculate(self):
+        positive_array = np.zeros((self.no_xy, self.no_xy, self.no_z_with_buffer))
+        negative_array = np.zeros((self.no_xy, self.no_xy, self.no_z_with_buffer))
+        positive_array_temp = np.zeros((self.no_xy, self.no_xy, self.no_z_with_buffer))
+        negative_array_temp = np.zeros((self.no_xy, self.no_xy, self.no_z_with_buffer))
+        no_recombined_charge_carriers = 0.0
+        no_initialised_charge_carriers = 0.0
+
+        """
+        Create an array which includes the number of track to be inserted for each time step
+        The tracks are distributed uniformly in time
+        """
+
+        positive_temp_entry = negative_temp_entry = recomb_temp = 0.0
+
+        f_steps_list = np.zeros(self.computation_time_steps)
+
+        szcz_pos = self.sz + self.cz * (self.cz + 1.0) / 2.0
+        szcz_neg = self.sz + self.cz * (self.cz - 1.0) / 2.0
+
+        sycy_pos = self.sy + self.cy * (self.cy + 1.0) / 2.0
+        sycy_neg = self.sy + self.cy * (self.cy - 1.0) / 2.0
+
+        sxcx_pos = self.sx + self.cx * (self.cx + 1.0) / 2.0
+        sxcx_neg = self.sx + self.cx * (self.cx - 1.0) / 2.0
+
+        cxyzsyz = (
+            1.0
+            - self.cx * self.cx
+            - self.cy * self.cy
+            - self.cz * self.cz
+            - 2.0 * (self.sx + self.sy + self.sz)
+        )
+
+        """
+        Start the simulation by evolving the distribution one step at a time
+        """
+
+        for time_step in range(self.computation_time_steps):
+
+            """
+            Refill the array with the electron density each time step
+            """
+            if self.should_simulate_beam_for_time_step(time_step):
+                no_initialised_step = self.simulate_beam(
+                    positive_array,
+                    negative_array,
+                )
+            else:
+                no_initialised_step = 0.0
+
+            no_initialised_charge_carriers += no_initialised_step
+
+            # calculate the new densities and store them in temporary arrays
+            for i in range(1, self.no_xy - 1):
+                for j in range(1, self.no_xy - 1):
+                    for k in range(1, self.no_z_with_buffer - 1):
+                        # using the Lax-Wendroff scheme
+                        positive_temp_entry = szcz_pos * positive_array[i, j, k - 1]
+                        positive_temp_entry += szcz_neg * positive_array[i, j, k + 1]
+
+                        positive_temp_entry += sycy_pos * positive_array[i, j - 1, k]
+                        positive_temp_entry += sycy_neg * positive_array[i, j + 1, k]
+
+                        positive_temp_entry += sxcx_pos * positive_array[i - 1, j, k]
+                        positive_temp_entry += sxcx_neg * positive_array[i + 1, j, k]
+
+                        positive_temp_entry += cxyzsyz * positive_array[i, j, k]
+
+                        # same for the negative charge carriers
+                        negative_temp_entry = szcz_pos * negative_array[i, j, k + 1]
+                        negative_temp_entry += szcz_neg * negative_array[i, j, k - 1]
+
+                        negative_temp_entry += sycy_pos * negative_array[i, j + 1, k]
+                        negative_temp_entry += sycy_neg * negative_array[i, j - 1, k]
+
+                        negative_temp_entry += sxcx_pos * negative_array[i + 1, j, k]
+                        negative_temp_entry += sxcx_neg * negative_array[i - 1, j, k]
+
+                        negative_temp_entry += cxyzsyz * negative_array[i, j, k]
+
+                        # the recombination part
+                        recomb_temp = (
+                            self.alpha
+                            * positive_array[i, j, k]
+                            * negative_array[i, j, k]
+                            * self.dt
+                        )
+                        positive_array_temp[i, j, k] = positive_temp_entry - recomb_temp
+                        negative_array_temp[i, j, k] = negative_temp_entry - recomb_temp
+
+                        if k > self.no_z_electrode and k < (
+                            self.no_z + self.no_z_electrode
+                        ):
+                            # sum over the recombination between the virtual electrodes
+                            no_recombined_charge_carriers += recomb_temp
+
+            f_steps_list[time_step] = (
+                no_initialised_charge_carriers - no_recombined_charge_carriers
+            ) / no_initialised_charge_carriers
+            np.copyto(positive_array, positive_array_temp)
+            np.copyto(negative_array, negative_array_temp)
+
+        return f_steps_list
diff --git a/electrons/python/pulsed_e_beam.py b/electrons/python/pulsed_e_beam.py
@@ -0,0 +1,34 @@
+from __future__ import division
+from generic_electron_solver import GenericElectronSolver
+
+
+class PulsedBeamPDEsolver(GenericElectronSolver):
+
+    def get_initialised_charge_carriers_after_beam(self):
+        initialised_charge_carriers = (
+            self.electron_density_per_cm3
+            * (self.no_xy - 2 * self.delta_border) ** 2
+            * self.no_z
+        )
+
+        return initialised_charge_carriers
+
+    def update_electron_density_arrays_after_beam(
+        self,
+        positive_array,
+        negative_array,
+    ):
+        positive_array[
+            self.delta_border : (self.no_xy - self.delta_border),
+            self.delta_border : (self.no_xy - self.delta_border),
+            self.no_z_electrode : (self.no_z + self.no_z_electrode),
+        ] += self.electron_density_per_cm3
+
+        negative_array[
+            self.delta_border : (self.no_xy - self.delta_border),
+            self.delta_border : (self.no_xy - self.delta_border),
+            self.no_z_electrode : (self.no_z + self.no_z_electrode),
+        ] += self.electron_density_per_cm3
+
+    def should_simulate_beam_for_time_step(self, time_step):
+        return time_step == 0