LPSE2D Bugs (#84)

* Factor of 2 in ky because of ny * intensity bug * updated config
ergodicio · Oct 25, 2024 · 6fcd5a8 · 6fcd5a8
1 parent d6108c3
commit 6fcd5a8
Show file tree

Hide file tree

Showing 4 changed files with 25 additions and 212 deletions.
diff --git a/adept/_lpse2d/core/epw.py b/adept/_lpse2d/core/epw.py
@@ -97,7 +97,6 @@ def tpd(self, t: float, y: Array, args: Dict) -> Array:
 
         tpd1 = E0[..., 1] * jnp.conj(ey)
         tpd1 = jnp.fft.ifft2(jnp.fft.fft2(tpd1) * self.low_pass_filter)
-        # tpd1 = E0_Ey
 
         divE_true = jnp.fft.ifft2(self.k_sq * jnp.fft.fft2(phi))
         E0_divE_k = jnp.fft.fft2(E0[..., 1] * jnp.conj(divE_true))
@@ -110,49 +109,6 @@ def tpd(self, t: float, y: Array, args: Dict) -> Array:
 
         return dphi
 
-    def calc_tpd1(self, t: float, y: Array, args: Dict) -> Array:
-        """
-        Calculates the first term of the two plasmon decay
-
-        Args:
-            t (float): time
-            y (Array): phi(x, y)
-            args (Dict): dictionary containing E0
-
-        Returns:
-            Array: dphi(x, y)
-        """
-        E0 = args["E0"]
-        phi = y
-
-        _, ey = self.calc_fields_from_phi(phi)
-
-        tpd1 = E0[..., 1] * jnp.conj(ey)
-        return self.tpd_const * tpd1
-
-    def calc_tpd2(self, t: float, y: Array, args: Dict) -> Array:
-        """
-        Calculates the second term of the two plasmon decay
-
-        Args:
-            t (float): time
-            y (Array): phi(x, y)
-            args (Dict): dictionary containing E0
-
-        Returns:
-            Array: dphi(x, y)
-
-        """
-        phi = y
-        E0 = args["E0"]
-
-        divE_true = jnp.fft.ifft2(self.k_sq * jnp.fft.fft2(phi))
-        E0_divE_k = jnp.fft.fft2(E0[..., 1] * jnp.conj(divE_true))
-
-        tpd2 = 1j * self.ky[None, :] * self.one_over_ksq * E0_divE_k
-        tpd2 = jnp.fft.ifft2(tpd2)
-        return self.tpd_const * tpd2
-
     def get_noise(self):
         random_amps = 1.0  # jax.random.uniform(self.amp_key, (self.nx, self.ny))
         random_phases = 2 * np.pi * jax.random.uniform(self.phase_key, (self.nx, self.ny))
@@ -175,8 +131,8 @@ def __call__(self, t: float, y: Dict[str, Array], args: Dict) -> Array:
         # density gradient
         if self.cfg["terms"]["epw"]["density_gradient"]:
             ex, ey = self.calc_fields_from_phi(phi)
-            ex *= jnp.exp(-1j * self.wp0 / 2.0 * (1 - background_density / self.envelope_density) * self.dt)
-            ey *= jnp.exp(-1j * self.wp0 / 2.0 * (1 - background_density / self.envelope_density) * self.dt)
+            ex *= jnp.exp(-1j * self.wp0 / 2.0 * (background_density / self.envelope_density - 1) * self.dt)
+            ey *= jnp.exp(-1j * self.wp0 / 2.0 * (background_density / self.envelope_density - 1) * self.dt)
             phi = self.calc_phi_from_fields(ex, ey)
 
         if self.cfg["terms"]["epw"]["source"]["noise"]:

diff --git a/adept/_lpse2d/core/integrator.py b/adept/_lpse2d/core/integrator.py
diff --git a/adept/_lpse2d/helpers.py b/adept/_lpse2d/helpers.py
@@ -34,9 +34,7 @@ def write_units(cfg: Dict) -> Dict:
     Z = cfg["units"]["ionization state"]
     A = cfg["units"]["atomic number"]
     lam0 = _Q(cfg["units"]["laser_wavelength"]).to("um").value
-    I0 = (
-        _Q(cfg["units"]["laser intensity"]).to("W/cm^2").value / 2
-    )  ## NB - the factor of 2 enables matching to the Short scaling
+    I0 = _Q(cfg["units"]["laser intensity"]).to("W/cm^2").value
     envelopeDensity = cfg["units"]["envelope density"]
 
     # Scaled constants
@@ -185,26 +183,16 @@ def get_derived_quantities(cfg: Dict) -> Dict:
     cfg_grid["ymin"] = _Q(cfg_grid["ymin"]).to("um").value
     cfg_grid["dx"] = _Q(cfg_grid["dx"]).to("um").value
 
-    cfg_grid["nx"] = int(cfg_grid["xmax"] / cfg_grid["dx"]) + 1
-    # cfg_grid["dx"] = cfg_grid["xmax"] / cfg_grid["nx"]
-    cfg_grid["dy"] = cfg_grid["dx"]  # cfg_grid["ymax"] / cfg_grid["ny"]
-    cfg_grid["ny"] = int(cfg_grid["ymax"] / cfg_grid["dy"]) + 1
+    cfg_grid["nx"] = int((cfg_grid["xmax"] - cfg_grid["xmin"]) / cfg_grid["dx"]) + 1
 
-    # midpt = (cfg_grid["xmax"] + cfg_grid["xmin"]) / 2
+    cfg_grid["dy"] = cfg_grid["dx"]
+    cfg_grid["ny"] = int((cfg_grid["ymax"] - cfg_grid["ymin"]) / cfg_grid["dy"]) + 1
 
-    # max_density = cfg["density"]["val at center"] + (cfg["grid"]["xmax"] - midpt) / L
-
-    # norm_n_max = np.abs(nmax / cfg["units"]["envelope density"] - 1)
-    # cfg_grid["dt"] = 0.05 * cfg_grid["dx"]
     cfg_grid["dt"] = _Q(cfg_grid["dt"]).to("ps").value
     cfg_grid["tmax"] = _Q(cfg_grid["tmax"]).to("ps").value
     cfg_grid["nt"] = int(cfg_grid["tmax"] / cfg_grid["dt"] + 1)
     cfg_grid["tmax"] = cfg_grid["dt"] * cfg_grid["nt"]
 
-    # if cfg_grid["nt"] > 1e6:
-    #     cfg_grid["max_steps"] = int(1e6)
-    #     print(r"Only running $10^6$ steps")
-    # else:
     cfg_grid["max_steps"] = cfg_grid["nt"] + 4
 
     # change driver parameters to the right units
@@ -308,7 +296,7 @@ def get_solver_quantities(cfg: Dict) -> Dict:
     )
 
     cfg_grid["zero_mask"] = (
-        np.where(cfg_grid["kx"] == 0, 0, 1)[:, None] * np.ones_like(cfg_grid["ky"])[None, :]
+        np.where(np.sqrt(cfg_grid["kx"][:, None] ** 2 + cfg_grid["ky"][None, :] ** 2) == 0, 0, 1)
         if cfg["terms"]["zero_mask"]
         else 1
     )

diff --git a/configs/envelope-2d/tpd.yaml b/configs/envelope-2d/tpd.yaml
@@ -1,7 +1,7 @@
 density:
   basis: linear
-  gradient scale length: 200.0um
-  max: 0.3
+  gradient scale length: 62.5um
+  max: 0.28
   min: 0.2
   noise:
     max: 1.0e-09
@@ -14,39 +14,41 @@ drivers:
     delta_omega_max: 0.015
     num_colors: 1
     envelope:
-      tw: 40ps
+      tw: 50ps
       tr: 0.25ps
-      tc: 20.5ps
+      tc: 20.05ps
       xr: 0.2um
       xw: 1000um
       xc: 50um
       yr: 0.2um
       yw: 1000um
       yc: 50um
+    params:
+      {}
 
 grid:
   boundary_abs_coeff: 1.0e4
-  boundary_width: 1um
-  low_pass_filter: 0.3
-  dt: 0.002ps
-  dx: 20nm
-  tmax: 25.0ps
+  boundary_width: 1.5um
+  low_pass_filter: 0.66
+  dt: 0.0005ps
+  dx: 40nm
+  tmax: 1.5ps
   tmin: 0.0ns
   ymax: 3um
   ymin: -3um
 mlflow:
   experiment: tpd
-  run: 1.5e14
+  run: test
 save:
   fields:
     t:
-      dt: 1ps
-      tmax: 25ps
+      dt: 0.1ps
+      tmax: 1.5ps
       tmin: 0ps
     x:
-      dx: 80nm
+      dx: 50nm
     y:
-      dy: 128nm
+      dy: 50nm
 solver: envelope-2d
 terms:
   epw:
@@ -62,11 +64,12 @@ terms:
       noise: true
       tpd: true
   zero_mask: true
+
 units:
   atomic number: 40
   envelope density: 0.25
   ionization state: 6
-  laser intensity: 2.0e+14W/cm^2
+  laser intensity: 1.2e+15W/cm^2
   laser_wavelength: 351nm
   reference electron temperature: 2000.0eV
   reference ion temperature: 1000eV