Add turbulence example

Author: FRidh

default.nix

@@ -16,6 +16,7 @@
 , pylint
 , pytest
 , pytestCheckHook
+, turbulence
 }:
 
 buildPythonPackage rec {
@@ -46,6 +47,7 @@ buildPythonPackage rec {
     pylint
     pytest
     pytestCheckHook
+    turbulence
   ];
 
   preCheck = ''

examples/basic_example.ipynb

@@ -96,7 +96,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = PSTD(maximum_frequency=maximum_frequency_target, pml=pml, medium=medium, cfl=0.05, size=[x, y])"
+    "model = PSTD(\n",
+    "    maximum_frequency=maximum_frequency_target,\n",
+    "    pml=pml,\n",
+    "    medium=medium,\n",
+    "    cfl=None,\n",
+    "    size=[x, y]\n",
+    ")"
    ]
   },
   {
@@ -221,10 +227,13 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
-    "model.run(seconds=0.060)"
+    "%%prun\n",
+    "model.run(seconds=0.06)"
    ]
   },
   {

examples/turbulence_example.ipynb

@@ -19,6 +19,7 @@
     "from pstd import PSTD, PML, Medium, Position2D, PointSource\n",
     "from pstd import PSTD\n",
     "from acoustics import Signal\n",
+    "from turbulence import Field2D, Gaussian2DTemp\n",
     "#import seaborn as sns\n",
     "%matplotlib inline"
    ]
@@ -38,14 +39,65 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "x = 30.0\n",
-    "y = 20.0\n",
+    "x = 50.0\n",
+    "y = 40.0\n",
     "z = 0.0\n",
     "\n",
-    "c = 343.2\n",
+    "c_0 = 343.2\n",
     "maximum_frequency_target = 200.0"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Turbulent field"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f_max = 500.0\n",
+    "\n",
+    "f_margin = 1.0\n",
+    "\n",
+    "# Amount of modes\n",
+    "max_mode_order = 100\n",
+    "# Maximum wavenumber\n",
+    "k_max = 10.0\n",
+    "\n",
+    "wavenumber_resolution = k_max / max_mode_order\n",
+    "\n",
+    "# We don't need it for the calculations but we do need it to create an instance.\n",
+    "spatial_resolution = c_0 / (2.0 * f_max * f_margin)\n",
+    "\n",
+    "spectrum = Gaussian2DTemp(\n",
+    "    max_mode_order=max_mode_order,\n",
+    "    wavenumber_resolution=wavenumber_resolution,\n",
+    "    mu_0=3e-2,\n",
+    "    a=1.1,\n",
+    "#     a=0.001,\n",
+    ")\n",
+    "\n",
+    "field = Field2D(\n",
+    "    x=x,\n",
+    "    y=y,\n",
+    "    z=y,\n",
+    "    spatial_resolution=spatial_resolution,\n",
+    "    spectrum=spectrum,\n",
+    ")\n",
+    "\n",
+    "mu = field.randomize().generate().mu\n",
+    "\n",
+    "print(\"Mu shape: {}\".format(mu.shape))\n",
+    "c = ( mu + 1.0 ) * c_0\n",
+    "\n",
+    "field.plot()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -93,7 +145,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = PSTD(maximum_frequency=maximum_frequency_target, pml=pml, medium=medium, cfl=0.05, size=[x, y])"
+    "model = PSTD(\n",
+    "    maximum_frequency=maximum_frequency_target,\n",
+    "    pml=pml,\n",
+    "    medium=medium,\n",
+    "    cfl=PSTD.maximum_cfl(medium.soundspeed)/2.,\n",
+    "    size=[x, y],\n",
+    "    spacing = spatial_resolution,\n",
+    ")"
    ]
   },
   {
@@ -109,7 +168,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "source_position = Position2D(x/4.0, y/2.0)\n",
+    "source_position = Position2D(x*2.0/5.0, y/2.0)\n",
     "\n",
     "source = model.add_object('source', 'PointSource', position=source_position, \n",
     "                         excitation='pulse', quantity='pressure', amplitude=0.1)"
@@ -128,7 +187,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "receiver_position = Position2D(x*3.0/4.0, y/2.0)\n",
+    "receiver_position = Position2D(x*3.0/5.0, y/2.0)\n",
     "\n",
     "receiver = model.add_object('receiver', 'Receiver', position=receiver_position, quantity='pressure')"
    ]
@@ -189,7 +248,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.run(steps=10)"
+    "model.run(seconds=0.01)"
    ]
   },
   {
@@ -221,14 +280,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.run(steps=500)"
+    "model.run(seconds=0.03)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "as you can see."
+    "With more steps we now see the effect of the turbulence."
    ]
   },
   {

flake.lock

@@ -3,8 +3,9 @@
 
   inputs.nixpkgs.url = "nixpkgs/nixpkgs-unstable";
   inputs.utils.url = "github:numtide/flake-utils";
+  inputs.turbulence.url = "github:FRidh/turbulence";
 
-  outputs = { self, nixpkgs, utils }: {
+  outputs = { self, nixpkgs, utils, turbulence }: {
     overlay = final: prev: {
       pythonPackagesOverrides = (prev.pythonPackagesOverrides or []) ++ [
         (self: super: {
@@ -23,7 +24,7 @@
   } // (utils.lib.eachSystem [ "x86_64-linux" ] (system: let
     pkgs = (import nixpkgs {
       inherit system;
-      overlays = [ self.overlay ];
+      overlays = [ self.overlay turbulence.overlay ];
     });
     python = pkgs.python3;
     pstd = python.pkgs.pstd;
@@ -42,4 +43,4 @@
       '';
     };
   }));
-}
+}

flake.nix

@@ -1220,7 +1220,7 @@ def __init__(
         maximum_frequency,
         medium=None,
         pml=None,
-        cfl=0.05,
+        cfl=None,
         spacing=None,
         axes=None,
         size=None,
@@ -1347,11 +1347,17 @@ def cfl(self):
         return self._cfl
 
     @cfl.setter
-    def cfl(self, x):
-        if not x <= np.mean(self.medium.soundspeed) / np.max(self.medium.soundspeed):
+    def cfl(self, x=None):
+        if x is None:
+            x = self.maximum_cfl(self.medium.soundspeed)
+        if not x <= self.maximum_cfl(self.medium.soundspeed):
             raise ValueError("CFL too high.")
         self._cfl = x
 
+    @staticmethod
+    def maximum_cfl(soundspeed):
+        return np.mean(soundspeed) / np.max(soundspeed)
+
     @property
     def constant_field(self):
         """Constant refractive-index field."""