use arrays for pml params

fdtd.py

@@ -156,8 +156,12 @@ def main():
         print('iteration {}: average {} iterations per sec'.format(t, (t+1)/(time.perf_counter()-start)))
         sys.stdout.flush()
 
+        if t % 20 == 0:
+            r = sum([(f * f * e).get().sum() for f, e in zip(sim.E, sim.eps)])
+            print('E sum', r)
+
         # Save field slices
-        if (t % 20 == 0 and (max_t - t < 1000 or t < 1000)) or t == max_t-1:
+        if (t % 20 == 0 and (max_t - t <= 1000 or t <= 2000)) or t == max_t-1:
             print('saving E-field')
             for j, f in enumerate(sim.E):
                 a = f.get()

fdtd/boundary.py

@@ -118,18 +118,12 @@ def cpml(direction: int,
         p0 = numpy.exp(-(sigma + alpha) * dt)
         p1 = sigma / (sigma + alpha) * (p0 - 1)
         return p0, p1
-    p0e, p1e = par(xe)
-    p0h, p1h = par(xh)
 
     vals = {'r': r,
             'u': uv[0],
             'v': uv[1],
             'np':   np,
             'th':   thickness,
-            'p0e': ', '.join((str(x) for x in p0e)),
-            'p1e': ', '.join((str(x) for x in p1e)),
-            'p0h': ', '.join((str(x) for x in p0h)),
-            'p1h': ', '.join((str(x) for x in p1h)),
             'se': '-+'[direction % 2],
             'sh': '+-'[direction % 2]}
 
@@ -149,24 +143,16 @@ if ( (s{r} - 1) > {r} && {r} > (s{r} - 1) - ({th} + 1) ) {{
         raise Exception('Bad polarity (=0)')
 
     code_e = """
-    // pml parameters:
-    const float p0[{th}] = {{ {p0e} }};
-    const float p1[{th}] = {{ {p1e} }};
-
-    Psi_{r}{np}_E{u}[ip] = p0[ir] * Psi_{r}{np}_E{u}[ip] + p1[ir] * (H{v}[i] - H{v}[i-di{r}]);
-    Psi_{r}{np}_E{v}[ip] = p0[ir] * Psi_{r}{np}_E{v}[ip] + p1[ir] * (H{u}[i] - H{u}[i-di{r}]);
+    Psi_{r}{np}_E{u}[ip] = p0e[ir] * Psi_{r}{np}_E{u}[ip] + p1e[ir] * (H{v}[i] - H{v}[i-di{r}]);
+    Psi_{r}{np}_E{v}[ip] = p0e[ir] * Psi_{r}{np}_E{v}[ip] + p1e[ir] * (H{u}[i] - H{u}[i-di{r}]);
 
     E{u}[i] {se}= dt / eps{u}[i] * Psi_{r}{np}_E{u}[ip];
     E{v}[i] {sh}= dt / eps{v}[i] * Psi_{r}{np}_E{v}[ip];
 }}
 """
     code_h = """
-    // pml parameters:
-    const float p0[{th}] = {{ {p0h} }};
-    const float p1[{th}] = {{ {p1h} }};
-
-    Psi_{r}{np}_H{u}[ip] = p0[ir] * Psi_{r}{np}_H{u}[ip] + p1[ir] * (E{v}[i+di{r}] - E{v}[i]);
-    Psi_{r}{np}_H{v}[ip] = p0[ir] * Psi_{r}{np}_H{v}[ip] + p1[ir] * (E{u}[i+di{r}] - E{u}[i]);
+    Psi_{r}{np}_H{u}[ip] = p0h[ir] * Psi_{r}{np}_H{u}[ip] + p1h[ir] * (E{v}[i+di{r}] - E{v}[i]);
+    Psi_{r}{np}_H{v}[ip] = p0h[ir] * Psi_{r}{np}_H{v}[ip] + p1h[ir] * (E{u}[i+di{r}] - E{u}[i]);
 
     H{u}[i] {sh}= dt * Psi_{r}{np}_H{u}[ip];
     H{v}[i] {se}= dt * Psi_{r}{np}_H{v}[ip];
@@ -180,6 +166,8 @@ if ( (s{r} - 1) > {r} && {r} > (s{r} - 1) - ({th} + 1) ) {{
                   'Psi_{r}{np}_E{v}'.format(**vals)],
         'psi_H': ['Psi_{r}{np}_H{u}'.format(**vals),
                   'Psi_{r}{np}_H{v}'.format(**vals)],
+        'pe': par(xe),
+        'ph': par(xh),
     }
 
     return pml_data

fdtd/simulation.py

@@ -171,14 +171,19 @@ class Simulation(object):
         dt_arg = [ctype + ' dt']
         arglist_E = [ctype + ' *' + psi for psi in psi_E_names]
         arglist_H = [ctype + ' *' + psi for psi in psi_H_names]
-        pml_E_args = ', '.join(E_args + H_args + dt_arg + eps_args + arglist_E)
-        pml_H_args = ', '.join(E_args + H_args + dt_arg + arglist_H)
+        pe_args = [ctype + ' *' + s for s in ('p0e', 'p1e')]
+        ph_args = [ctype + ' *' + s for s in ('p0h', 'p1h')]
+        pml_E_args = ', '.join(E_args + H_args + dt_arg + eps_args + arglist_E + pe_args)
+        pml_H_args = ', '.join(E_args + H_args + dt_arg + arglist_H + ph_args)
 
         pml_E = ElementwiseKernel(self.context, arguments=pml_E_args, operation=pml_E_source)
         pml_H = ElementwiseKernel(self.context, arguments=pml_H_args, operation=pml_H_source)
 
-        self.cpml_E = lambda e: pml_E(*self.E, *self.H, self.dt, *self.eps, *psi_E, wait_for=e)
-        self.cpml_H = lambda e: pml_H(*self.E, *self.H, self.dt, *psi_H, wait_for=e)
+        pe = [pyopencl.array.to_device(self.queue, p) for p in pml_data['pe']]
+        ph = [pyopencl.array.to_device(self.queue, p) for p in pml_data['ph']]
+
+        self.cpml_E = lambda e: pml_E(*self.E, *self.H, self.dt, *self.eps, *psi_E, *pe, wait_for=e)
+        self.cpml_H = lambda e: pml_H(*self.E, *self.H, self.dt, *psi_H, *ph, wait_for=e)
         self.cmpl_psi_E = {k: v for k, v in zip(psi_E_names, psi_E)}
         self.cmpl_psi_H = {k: v for k, v in zip(psi_H_names, psi_H)}