Merge pull request #5 from uwoseis/windtunnel

Merge in current changes from branch 'windtunnel'

anemoi/Tests/test_Analytical.py

@@ -1,7 +1,7 @@
 
 import unittest
 import numpy as np
-from anemoi import AnalyticalHelmholtz
+from anemoi import AnalyticalHelmholtz, FakeSource
 
 class TestAnalyticalHelmholtz(unittest.TestCase):
 
@@ -19,9 +19,11 @@ def test_cleanExecution(self):
             'nz':       nz,     # count
             'freq':     2e2,    # Hz
         }
-
+
+        sloc = np.array([nx/2, nz/2]).reshape((1,2))
+
         Green = AnalyticalHelmholtz(systemConfig)
-        u = Green(nx/2, nz/2)
+        u = Green(sloc)
 
 if __name__ == '__main__':
     unittest.main()

anemoi/Tests/test_Eurus.py

@@ -31,7 +31,8 @@ def test_cleanExecution(self):
         Ainv = Eurus(systemConfig)
         src = StackedSimpleSource(systemConfig)
 
-        q = src(nx/2, nz/2)
+        sloc = np.array([nx/2, nz/2]).reshape((1,2))
+        q = src(sloc)
         u = Ainv*q
 
     def test_compareAnalytical_Isotropic(self):
@@ -68,14 +69,15 @@ def test_compareAnalytical_Isotropic(self):
 
         xs = 25
         zs = 25
+        sloc = np.array([xs, zs]).reshape((1,2))
 
         Ainv = Eurus(systemConfig)
         src = StackedSimpleSource(systemConfig)
-        q = src(xs, zs)
+        q = src(sloc)
         uMZ = Ainv*q
 
         AH = AnalyticalHelmholtz(systemConfig)
-        uAH = AH(xs, zs)
+        uAH = AH(sloc)
 
         nx = systemConfig['nx']
         nz = systemConfig['nz']
@@ -124,14 +126,15 @@ def test_compareAnalytical_Elliptical(self):
 
         xs = 25
         zs = 25
+        sloc = np.array([xs, zs]).reshape((1,2))
 
         Ainv = Eurus(systemConfig)
         src = StackedSimpleSource(systemConfig)
-        q = src(xs, zs)
+        q = src(sloc)
         uMZ = Ainv*q
 
         AH = AnalyticalHelmholtz(systemConfig)
-        uAH = AH(xs, zs)
+        uAH = AH(sloc)
 
         nx = systemConfig['nx']
         nz = systemConfig['nz']

anemoi/Tests/test_MiniZephyr.py

@@ -23,11 +23,12 @@ def test_cleanExecution(self):
 
         xs = 50
         zs = 100
+        sloc = np.array([xs, zs]).reshape((1,2))
 
         Ainv = MiniZephyr(systemConfig)
         src = SimpleSource(systemConfig)
 
-        q = src(xs, zs)
+        q = src(sloc)
         u = Ainv*q
 
     def test_cleanExecution25D(self):
@@ -44,11 +45,12 @@ def test_cleanExecution25D(self):
 
         xs = 50
         zs = 100
+        sloc = np.array([xs, zs]).reshape((1,2))
 
         Ainv = MiniZephyr25D(systemConfig)
         src = SimpleSource(systemConfig)
 
-        q = src(xs, zs)
+        q = src(sloc)
         u = Ainv*q
 
     def test_cleanExecution25DParallel(self):
@@ -65,11 +67,12 @@ def test_cleanExecution25DParallel(self):
 
         xs = 50
         zs = 100
+        sloc = np.array([xs, zs]).reshape((1,2))
 
         Ainv = MiniZephyr25D(systemConfig)
         src = SimpleSource(systemConfig)
 
-        q = src(xs, zs)
+        q = src(sloc)
         u = Ainv*q
 
     def test_compareAnalytical(self):
@@ -84,14 +87,15 @@ def test_compareAnalytical(self):
 
         xs = 25 
         zs = 25
+        sloc = np.array([xs, zs]).reshape((1,2))
 
         Ainv = MiniZephyr(systemConfig)
         src = SimpleSource(systemConfig)
-        q = src(xs, zs)
+        q = src(sloc)
         uMZ = Ainv*q
 
         AH = AnalyticalHelmholtz(systemConfig)
-        uAH = AH(xs, zs)
+        uAH = AH(sloc)
 
         nx = systemConfig['nx']
         nz = systemConfig['nz']
@@ -120,14 +124,15 @@ def test_compareAnalytical25D(self):
 
         xs = 25 
         zs = 25
+        sloc = np.array([xs, zs]).reshape((1,2))
 
         Ainv = MiniZephyr25D(systemConfig)
         src = SimpleSource(systemConfig)
-        q = src(xs, zs)
+        q = src(sloc)
         uMZ = Ainv*q
 
         AH = AnalyticalHelmholtz(systemConfig)
-        uAH = AH(xs, zs)
+        uAH = AH(sloc)
 
         nx = systemConfig['nx']
         nz = systemConfig['nz']

anemoi/Tests/test_Sources.py

@@ -0,0 +1,68 @@
+import unittest
+import numpy as np
+from anemoi import SimpleSource, StackedSimpleSource, SparseKaiserSource, KaiserSource
+
+class TestSources(unittest.TestCase):
+
+    @staticmethod
+    def elementNorm(arr):
+        return np.sqrt((arr.conj()*arr).sum()) / arr.size
+
+    def setUp(self):
+        pass
+
+    def test_cleanExecution(self):
+
+        systemConfig = {
+            'nx':       100,
+            'nz':       100,
+        }
+
+        loc = np.array([[50.,50.],[25.,25.], [80.,80.], [25.,80.]])
+
+        ss  = SimpleSource(systemConfig)
+        sss = SimpleSource(systemConfig)
+        sks = SparseKaiserSource(systemConfig)
+        ks  = KaiserSource(systemConfig)
+
+        qss     = ss(loc)
+        qsss    = sss(loc)
+        qsks    = sks(loc)
+        qks     = ks(loc)
+
+    def test_KaiserSource(self):
+
+        systemConfig = {
+            'nx':       100,
+            'nz':       100,
+        }
+
+        loc = np.array([[50.,50.],[25.,25.], [80.,80.], [25.,80.]])
+
+        sks = SparseKaiserSource(systemConfig)
+        ks  = KaiserSource(systemConfig)
+
+        qsks    = sks(loc)
+        qks     = ks(loc)
+
+        self.assertTrue(self.elementNorm(qsks.toarray() - qks) == 0.)
+
+    def test_KaiserSourceSimpleCase(self):
+
+        systemConfig = {
+            'nx':       100,
+            'nz':       100,
+            'dx':       1.,  # NB: This will not generally be true,
+            'dz':       1.,  # but the scaleTerm in KS is 1/(dx*dz).
+        }
+
+        loc = np.array([[50.,50.],[25.,25.], [80.,80.], [25.,80.]])
+
+        ss  = SimpleSource(systemConfig)
+        ks  = KaiserSource(systemConfig)
+
+        qss = ss(loc)
+        qks = ks(loc)
+
+        self.assertTrue(self.elementNorm(qks - qss) < 1e-10)
+

anemoi/__init__.py

@@ -1,4 +1,7 @@
+
 from .analytical import AnalyticalHelmholtz
+from .discretization import BaseDiscretization, DiscretizationWrapper, MultiFreq
 from .minizephyr import MiniZephyr, MiniZephyr25D
 from .eurus import Eurus
-from .source import SimpleSource, StackedSimpleSource
+from .source import BaseSource, SimpleSource, StackedSimpleSource, KaiserSource, SparseKaiserSource, FakeSource
+from .meta import BaseModelDependent, BaseSCCache, AttributeMapper, SCFilter

anemoi/analytical.py

@@ -41,7 +41,10 @@ def Green3D(self, r):
         # Correct: (1./(4*np.pi*r)) * np.exp(-1j*self.k*r)
         return self.scaleterm * (1./(4*np.pi*r)) * np.exp(1j*self.k*r)
 
-    def __call__(self, x, z):
+    def __call__(self, q):
+
+        x = q[0,0]
+        z = q[0,-1]
 
         dx = self._x - x
         dz = self._z - z
@@ -52,3 +55,7 @@ def __call__(self, x, z):
         stretch = np.sqrt(self.stretch * np.cos(strangle)**2 + np.sin(strangle)**2)
 
         return np.nan_to_num(self.Green(dist * stretch)).ravel()
+
+    def __mul__(self, q):
+
+        return self(q)

anemoi/discretization.py

@@ -1,7 +1,17 @@
 
-from .meta import AttributeMapper, BaseModelDependent
-from .solver import DirectSolver
+import copy
 import numpy as np
+from .meta import AttributeMapper, BaseSCCache, BaseModelDependent
+from .solver import DirectSolver
+
+try:
+    from multiprocessing import Pool, Process
+except ImportError:
+    PARALLEL = False
+else:
+    PARALLEL = True
+
+PARTASK_TIMEOUT = 60
 
 class BaseDiscretization(BaseModelDependent):
 
@@ -42,3 +52,78 @@ def __mul__(self, rhs):
 
     def __call__(self, value):
         return self*value
+
+class DiscretizationWrapper(BaseSCCache):
+
+    initMap = {
+    #   Argument        Required    Rename as ...   Store as type
+        'disc':         (True,      None,           None),
+        'scaleTerm':    (False,     '_scaleTerm',   np.complex128),
+    }
+
+    @property
+    def scaleTerm(self):
+        return getattr(self, '_scaleTerm', 1.)
+
+    @property
+    def _spConfigs(self):
+
+        def duplicateUpdate(spu):
+            nsc = copy.copy(self.systemConfig)
+            nsc.update(spu)
+            return nsc
+
+        return (duplicateUpdate(spu) for spu in self.spUpdates)
+
+    @property
+    def subProblems(self):
+        if getattr(self, '_subProblems', None) is None:
+
+            self._subProblems = map(self.disc, self._spConfigs)
+        return self._subProblems
+
+    @property
+    def spUpdates(self):
+        raise NotImplementedError
+
+    def __mul__(self, rhs):
+        raise NotImplementedError
+
+
+class MultiFreq(DiscretizationWrapper):
+
+    initMap = {
+    #   Argument        Required    Rename as ...   Store as type
+        'disc':         (True,      '_disc',        None),
+        'freqs':        (True,      None,           list),
+        'parallel':     (False,     '_parallel',    bool),
+    }
+
+    maskKeys = ['disc', 'freqs', 'parallel']
+
+    @property
+    def parallel(self):
+        return PARALLEL and getattr(self, '_parallel', True)
+
+    @property
+    def spUpdates(self):
+        return [{'freq': freq} for freq in self.freqs]
+
+    @property
+    def disc(self):
+        return self._disc
+
+    def __mul__(self, rhs):
+
+        if self.parallel:
+            pool = Pool()
+            plist = []
+            for sp in self.subProblems:
+                p = pool.apply_async(sp, (rhs,))
+                plist.append(p)
+
+            u = (self.scaleTerm*p.get(PARTASK_TIMEOUT) for p in plist)
+        else:
+            u = (self.scaleTerm*(sp*rhs) for sp in self.subProblems)
+
+        return list(u) # TODO: Maybe we *do* want to return a generator here?