Compare commits
1 Commits
Author | SHA1 | Date | |
---|---|---|---|
727105f8fe |
6
fdtd.py
6
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()
|
||||
|
@ -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
|
||||
|
@ -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)}
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user