Compare commits

...

1 Commits

Author SHA1 Message Date
jan
727105f8fe use arrays for pml params 2016-09-01 14:39:44 -07:00
3 changed files with 20 additions and 23 deletions

View File

@ -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()

View File

@ -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

View File

@ -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)}