Fix model stop on speed up and test

Author: Sosnowsky

blobmodel/model.py

@@ -339,9 +339,9 @@ def _compute_start_stop(self, blob: Blob, speed_up: bool, error: float):
         # ignores t_drain when calculating stop time
         stop = np.minimum(
             self._geometry.t.size,
-            start
-            + int(
-                (
+            int(
+                blob.t_init
+                + (
                     -blob.width_prop * np.log(error * np.sqrt(np.pi))
                     + self._geometry.Lx
                     - blob.pos_x

examples/compare_to_analytical_sol.py

@@ -30,7 +30,7 @@
 )
 
 ds = tmp.make_realization(speed_up=True, error=1e-10)
-x = np.linspace(0, 10, 100)
+x = np.linspace(0, 10, 20)
 t_p = 1
 t_w = 1 / 10
 amp = 1

tests/test_analytical.py

@@ -7,51 +7,46 @@
 )
 import xarray as xr
 import numpy as np
+import matplotlib.pyplot as plt
 
 
 # use DefaultBlobFactory to define distribution functions fo random variables
 bf = DefaultBlobFactory(A_dist=DistributionEnum.deg, vy_dist=DistributionEnum.zeros)
 
-t_loss = 2.0
-
+t_drain = 2
+Nx, Lx = 10, 10
 tmp = Model(
-    Nx=100,
+    Nx=Nx,
     Ny=1,
-    Lx=10,
+    Lx=Lx,
     Ly=0,
-    dt=0.1,
+    dt=1,
     T=1000,
-    t_drain=t_loss,
-    blob_shape=BlobShapeImpl(
-        BlobShapeEnum.exp, BlobShapeEnum.gaussian
-    ),  # Analytical form only applies to exp shape
+    blob_shape=BlobShapeImpl(BlobShapeEnum.exp, BlobShapeEnum.gaussian),
+    t_drain=t_drain,
     periodic_y=False,
     num_blobs=10000,
     blob_factory=bf,
-    t_init=10,
     one_dimensional=True,
 )
 
-tmp.make_realization(file_name="test_analytical.nc", speed_up=True, error=1e-2)
-
 
 def test_convergence_to_analytical_solution():
     """
     Checks that the mean value of a one-dimensional realization with constant velocities agrees with the
     analytical derived results. See O. E. Garcia, et al.; Phys. Plasmas 1 May 2016; 23 (5): 052308. https://doi.org/10.1063/1.4951016
     """
-    ds = xr.open_dataset("test_analytical.nc")
+    ds = tmp.make_realization(speed_up=True, error=1e-10)
     model_profile = ds.n.isel(y=0).mean(dim="t")
 
-    x = np.linspace(0, 10, 100)
+    x = np.arange(0, Lx, Lx / Nx)
     t_p = 1  # vx/ell
     t_w = 1 / 10
     amp = 1
     v_p = 1.0
-    t_d = t_loss * t_p / (t_loss + t_p)
-
-    analytical_profile = t_d / t_w * amp * np.exp(-x / (v_p * t_loss))
+    t_d = t_drain * t_p / (t_drain + t_p)
 
+    analytical_profile = t_d / t_w * amp * np.exp(-x / (v_p * t_drain))
     error = np.mean(abs(model_profile.values - analytical_profile))
 
     assert error < 0.1, "Numerical error too big"