Merge pull request #33 from uit-cosmo/bug_vy=0

bugfix for periodic_y = True and vy = 0

blobmodel/blobs.py

@@ -65,20 +65,21 @@ def discretize_blob(
             origin=(self.pos_x, self.pos_y), x=x, y=y, angle=-self.theta
         )
         if periodic_y:
-            __x_border = (Ly - self.pos_y) / np.sin(self.theta)
+            if np.sin(self.theta) == 0:
+                __x_border = Ly - self.pos_y
+                __adjusted_Ly = Ly
+            else:
+                __x_border = (Ly - self.pos_y) / np.sin(self.theta)
+                __adjusted_Ly = Ly / np.sin(self.theta)
             if type(t) == int or type(t) == float:
                 # t has dimensionality = 0, used for testing
                 __number_of_y_propagations = (
-                    self.__prop_dir_blob_position(t)
-                    + Ly / np.sin(self.theta)
-                    - __x_border
-                ) // (Ly / np.sin(self.theta))
+                    self.__prop_dir_blob_position(t) + __adjusted_Ly - __x_border
+                ) // __adjusted_Ly
             else:
                 __number_of_y_propagations = (
-                    self.__prop_dir_blob_position(t)[0, 0]
-                    + Ly / np.sin(self.theta)
-                    - __x_border
-                ) // (Ly / np.sin(self.theta))
+                    self.__prop_dir_blob_position(t)[0, 0] + __adjusted_Ly - __x_border
+                ) // __adjusted_Ly
             __blob_values = (
                 self.__single_blob(
                     x_perp, y_perp, t, Ly, periodic_y, __number_of_y_propagations

blobmodel/model.py

@@ -111,15 +111,14 @@ def make_realization(
             # can also be used for gaussian pulses since they converge faster than exponential pulses
             if speed_up:
                 start = int(b.t_init / self.__geometry.dt)
-                try:
+                if b.v_x == 0:
+                    stop = self.__geometry.t.size
+                else:
                     # ignores t_drain when calculating stop time
                     stop = start + int(
                         (-np.log(error * np.sqrt(np.pi)) + self.__geometry.Lx - b.pos_x)
                         / (b.v_x * self.__geometry.dt)
                     )
-                except:
-                    print("Warning occurs due to v_x == 0")
-                    stop = self.__geometry.t.size
                 output[:, :, start:stop] += b.discretize_blob(
                     x=self.__geometry.x_matrix[:, :, start:stop],
                     y=self.__geometry.y_matrix[:, :, start:stop],

tests/test_vx_or_vy=0.py

@@ -0,0 +1,125 @@
+from blobmodel import Model, BlobFactory, Blob
+import numpy as np
+
+
+class CustomBlobFactoryVy0(BlobFactory):
+    def __init__(self) -> None:
+        pass
+
+    def sample_blobs(
+        self, Ly: float, T: float, num_blobs: int, blob_shape: str, t_drain: float
+    ) -> list[Blob]:
+
+        # set custom parameter distributions
+        __amp = np.ones(num_blobs)
+        __width = np.ones(num_blobs)
+        __vx = np.ones(num_blobs)
+        __vy = np.zeros(num_blobs)
+
+        __posx = np.zeros(num_blobs)
+        __posy = np.ones(num_blobs) * Ly / 2
+        __t_init = np.ones(num_blobs) * 0
+
+        # this block must remain the same
+        __blobs = []
+        for i in range(num_blobs):
+            __blobs.append(
+                Blob(
+                    id=i,
+                    blob_shape=blob_shape,
+                    amplitude=__amp[i],
+                    width_prop=__width[i],
+                    width_perp=__width[i],
+                    v_x=__vx[i],
+                    v_y=__vy[i],
+                    pos_x=__posx[i],
+                    pos_y=__posy[i],
+                    t_init=__t_init[i],
+                    t_drain=t_drain,
+                )
+            )
+        return __blobs
+
+
+class CustomBlobFactoryVx0(BlobFactory):
+    def __init__(self) -> None:
+        pass
+
+    def sample_blobs(
+        self, Ly: float, T: float, num_blobs: int, blob_shape: str, t_drain: float
+    ) -> list[Blob]:
+
+        # set custom parameter distributions
+        __amp = np.ones(num_blobs)
+        __width = np.ones(num_blobs)
+        __vx = np.zeros(num_blobs)
+        __vy = np.ones(num_blobs)
+
+        __posx = np.zeros(num_blobs)
+        __posy = np.ones(num_blobs) * Ly / 2
+        __t_init = np.ones(num_blobs) * 0
+
+        # this block must remain the same
+        __blobs = []
+        for i in range(num_blobs):
+            __blobs.append(
+                Blob(
+                    id=i,
+                    blob_shape=blob_shape,
+                    amplitude=__amp[i],
+                    width_prop=__width[i],
+                    width_perp=__width[i],
+                    v_x=__vx[i],
+                    v_y=__vy[i],
+                    pos_x=__posx[i],
+                    pos_y=__posy[i],
+                    t_init=__t_init[i],
+                    t_drain=t_drain,
+                )
+            )
+        return __blobs
+
+
+bf_vy_0 = CustomBlobFactoryVy0()
+
+bm_vy_0 = Model(
+    Nx=10,
+    Ny=10,
+    Lx=10,
+    Ly=10,
+    dt=1,
+    T=1,
+    periodic_y=True,
+    blob_shape="exp",
+    num_blobs=1,
+    blob_factory=bf_vy_0,
+    t_drain=1e10,
+)
+
+bf_vx_0 = CustomBlobFactoryVx0()
+
+bm_vx_0 = Model(
+    Nx=10,
+    Ny=10,
+    Lx=10,
+    Ly=10,
+    dt=1,
+    T=1,
+    periodic_y=True,
+    blob_shape="exp",
+    num_blobs=1,
+    blob_factory=bf_vx_0,
+    t_drain=1e10,
+)
+
+
+def test_vy_0():
+    assert bm_vy_0.make_realization(speed_up=True, error=1e-2)
+
+
+def test_vx_0():
+    assert bm_vx_0.make_realization(speed_up=True, error=1e-2)
+
+
+test_vx_0()
+test_vy_0()