From 727105f8fe114f62233820203b5c6eadbf1638ef Mon Sep 17 00:00:00 2001 From: jan Date: Thu, 1 Sep 2016 14:39:44 -0700 Subject: [PATCH] use arrays for pml params --- fdtd.py | 6 +++++- fdtd/boundary.py | 24 ++++++------------------ fdtd/simulation.py | 13 +++++++++---- 3 files changed, 20 insertions(+), 23 deletions(-) diff --git a/fdtd.py b/fdtd.py index 0ce560e..63d2337 100644 --- a/fdtd.py +++ b/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() diff --git a/fdtd/boundary.py b/fdtd/boundary.py index bb8dbdf..6b50f9f 100644 --- a/fdtd/boundary.py +++ b/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 diff --git a/fdtd/simulation.py b/fdtd/simulation.py index 227dcb0..6c63273 100644 --- a/fdtd/simulation.py +++ b/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)}