Merge pull request #72 from kaelyndunnell/gaussian_distribution_fix

Gaussian distribution fix

src/hisp/festim_models/mb_model.py

@@ -48,10 +48,11 @@ def make_W_mb_model(
 
     vertices = np.concatenate(  # 1D mesh with extra refinement
         [
-            np.linspace(0, 30e-9, num=200),
-            np.linspace(30e-9, 3e-6, num=300),
-            np.linspace(3e-6, 30e-6, num=200),
-            np.linspace(30e-6, L, num=200),
+            np.linspace(0, 30e-9, num=300),
+            np.linspace(30e-9, 3e-6, num=400),
+            np.linspace(3e-6, 30e-6, num=400),
+            np.linspace(30e-6, 1e-4, num=400),
+            np.linspace(1e-4, L, num=300),
         ]
     )
     my_model.mesh = F.Mesh1D(vertices)
@@ -274,7 +275,7 @@ def make_W_mb_model(
 
     ############# Settings #############
     my_model.settings = F.Settings(
-        atol=1e-10,
+        atol=1e10,
         rtol=1e-10,
         max_iterations=1000,
         final_time=final_time,
@@ -317,8 +318,9 @@ def make_B_mb_model(
 
     vertices = np.concatenate(  # 1D mesh with extra refinement
         [
-            np.linspace(0, 30e-9, num=200),
-            np.linspace(30e-9, L, num=200),
+            np.linspace(0, 5e-9, num=200),
+            np.linspace(5e-9, 1e-8, num=500),
+            np.linspace(1e-8, L, num=300),
         ]
     )
     my_model.mesh = F.Mesh1D(vertices)
@@ -538,8 +540,8 @@ def make_B_mb_model(
 
     ############# Settings #############
     my_model.settings = F.Settings(
-        atol=1e-15,
-        rtol=1e-15,
+        atol=1e10,
+        rtol=1e-8,
         max_iterations=1000,
         final_time=final_time,
     )
@@ -743,8 +745,8 @@ def make_DFW_mb_model(
 
     ############# Settings #############
     my_model.settings = F.Settings(
-        atol=1e-15,
-        rtol=1e-15,
+        atol=1e10,
+        rtol=1e-10,
         max_iterations=1000,
         final_time=final_time,
     )

src/hisp/helpers.py

@@ -40,7 +40,20 @@ def update(self, t: float):
 def gaussian_distribution(
     x: npt.NDArray, mean: float, width: float, mod=ufl
 ) -> ufl.core.expr.Expr:
-    return mod.exp(-((x[0] - mean) ** 2) / (2 * width**2))
+    """Generates a gaussian distribution for particle sources.
+
+    Args:
+        x (npt.NDArray): x values along the length of given bin.
+        mean (float): Mean of the distribution.
+        width (float): Width of the gaussian distribution.
+        mod (_type_, optional): Module used to express gaussian distribution. Defaults to ufl.
+
+    Returns:
+        ufl.core.expr.Expr: Gaussian distribution with area 1.  
+    """
+    return mod.exp(-((x[0] - mean) ** 2) / (2 * width**2)) / (
+        np.sqrt(2 * np.pi * width**2)
+    )
 
 
 def periodic_step_function(x, period_on, period_total, value, value_off=0.0):

src/hisp/plamsa_data_handling/main.py

@@ -223,15 +223,17 @@ def get_heat(self, pulse: Pulse, bin: SubBin | DivBin, t_rel: float) -> float:
             raise ValueError(f"Invalid pulse type {pulse.pulse_type}")
 
         if pulse.pulse_type == "FP":
-            heat_total = data["heat_total"][bin_index]
+            photon_heat_radiation = 0.11e6  # W/m2
+            heat_total = data["heat_total"][bin_index]+ photon_heat_radiation 
             heat_ion = data["heat_ion"][bin_index]
             if isinstance(bin, SubBin):
                 heat_val = heat_total - heat_ion * (1 - bin.wetted_frac)
             else:
                 heat_val = heat_total
         elif pulse.pulse_type == "RISP":
             if isinstance(bin, SubBin):
-                heat_total = data["heat_total"]
+                photon_radiation_heat = 0.11e6  # W/m2
+                heat_total = data["heat_total"]+photon_radiation_heat 
                 heat_ion = data["heat_ion"]
                 heat_val = heat_total - heat_ion * (1 - bin.wetted_frac)
             else:

test/test_helpers.py

@@ -0,0 +1,13 @@
+from hisp.helpers import gaussian_distribution
+
+import numpy as np
+
+
+def test_gaussian_distr_area():
+    """Generates a gaussian distribution and checks if the area under the curve is 1."""
+    x = np.linspace(-10, 10, 1000)
+    y = gaussian_distribution([x], mean=0, width=1, mod=np)
+    area = np.trapz(y, x)
+    assert np.isclose(
+        area, 1, atol=1e-2
+    ), f"Area under Gaussian distribution is {area} instead of 1"