forked from jan/opencl_fdtd
Rewrite, with the following features:
- Move to jinja2 templates for the opencl code - Combine PML code into the E, H updates for speed - Add Poynting vector calculation code, including precalculation during H update - Use arrays for PML parameters (p0, p1) - Switch to linearized, C-ordered fields (~50% performance boost??) - Added jinja2 and fdfd_tools dependenciesimage
jan
8 years ago
committed by
Jan Petykiewicz
parent
cd72219d0b
commit
d34c478f1d
@ -1,72 +0,0 @@
|
||||
"""
|
||||
Basic code snippets for opencl FDTD
|
||||
"""
|
||||
|
||||
from typing import List
|
||||
import numpy
|
||||
|
||||
|
||||
def shape_source(shape: List[int] or numpy.ndarray) -> str:
|
||||
"""
|
||||
Defines sx, sy, sz C constants specifying the shape of the grid in each of the 3 dimensions.
|
||||
|
||||
:param shape: [sx, sy, sz] values.
|
||||
:return: String containing C source.
|
||||
"""
|
||||
sxyz = """
|
||||
// Field sizes
|
||||
const int sx = {shape[0]};
|
||||
const int sy = {shape[1]};
|
||||
const int sz = {shape[2]};
|
||||
""".format(shape=shape)
|
||||
return sxyz
|
||||
|
||||
# Defines dix, diy, diz constants used for stepping in the x, y, z directions in a linear array
|
||||
# (ie, given Ex[i] referring to position (x, y, z), Ex[i+diy] will refer to position (x, y+1, z))
|
||||
dixyz_source = """
|
||||
// Convert offset in field xyz to linear index offset
|
||||
const int dix = sz * sy;
|
||||
const int diy = sz;
|
||||
const int diz = 1;
|
||||
"""
|
||||
|
||||
# Given a linear index i and shape sx, sy, sz, defines x, y, and z as the 3D indices of the current element (i).
|
||||
xyz_source = """
|
||||
// Convert linear index to field index (xyz)
|
||||
const int x = i / (sz * sy);
|
||||
const int y = (i - x * sz * sy) / sz;
|
||||
const int z = (i - y * sz - x * sz * sy);
|
||||
"""
|
||||
|
||||
# Source code for updating the E field; maxes use of dixyz_source.
|
||||
maxwell_E_source = """
|
||||
// E update equations
|
||||
Ex[i] += dt / epsx[i] * ((Hz[i] - Hz[i-diy]) - (Hy[i] - Hy[i-diz]));
|
||||
Ey[i] += dt / epsy[i] * ((Hx[i] - Hx[i-diz]) - (Hz[i] - Hz[i-dix]));
|
||||
Ez[i] += dt / epsz[i] * ((Hy[i] - Hy[i-dix]) - (Hx[i] - Hx[i-diy]));
|
||||
"""
|
||||
|
||||
# Source code for updating the H field; maxes use of dixyz_source and assumes mu=0
|
||||
maxwell_H_source = """
|
||||
// H update equations
|
||||
Hx[i] -= dt * ((Ez[i+diy] - Ez[i]) - (Ey[i+diz] - Ey[i]));
|
||||
Hy[i] -= dt * ((Ex[i+diz] - Ex[i]) - (Ez[i+dix] - Ez[i]));
|
||||
Hz[i] -= dt * ((Ey[i+dix] - Ey[i]) - (Ex[i+diy] - Ex[i]));
|
||||
"""
|
||||
|
||||
|
||||
def type_to_C(float_type: numpy.float32 or numpy.float64) -> str:
|
||||
"""
|
||||
Returns a string corresponding to the C equivalent of a numpy type.
|
||||
Only works for float32 and float64.
|
||||
|
||||
:param float_type: numpy.float32 or numpy.float64
|
||||
:return: string containing the corresponding C type (eg. 'double')
|
||||
"""
|
||||
if float_type == numpy.float32:
|
||||
arg_type = 'float'
|
||||
elif float_type == numpy.float64:
|
||||
arg_type = 'double'
|
||||
else:
|
||||
raise Exception('Unsupported type')
|
||||
return arg_type
|
||||
@ -1,185 +0,0 @@
|
||||
from typing import List, Dict
|
||||
import numpy
|
||||
|
||||
|
||||
def conductor(direction: int,
|
||||
polarity: int,
|
||||
) -> List[str]:
|
||||
"""
|
||||
Create source code for conducting boundary conditions.
|
||||
|
||||
:param direction: integer in range(3), corresponding to x,y,z.
|
||||
:param polarity: -1 or 1, specifying eg. a -x or +x boundary.
|
||||
:return: [E_source, H_source] source code for E and H boundary update steps.
|
||||
"""
|
||||
if direction not in range(3):
|
||||
raise Exception('Invalid direction: {}'.format(direction))
|
||||
|
||||
if polarity not in (-1, 1):
|
||||
raise Exception('Invalid polarity: {}'.format(polarity))
|
||||
|
||||
r = 'xyz'[direction]
|
||||
uv = 'xyz'.replace(r, '')
|
||||
|
||||
if polarity < 0:
|
||||
bc_e = """
|
||||
if ({r} == 0) {{
|
||||
E{r}[i] = 0;
|
||||
E{u}[i] = E{u}[i+di{r}];
|
||||
E{v}[i] = E{v}[i+di{r}];
|
||||
}}
|
||||
"""
|
||||
bc_h = """
|
||||
if ({r} == 0) {{
|
||||
H{r}[i] = H{r}[i+di{r}];
|
||||
H{u}[i] = 0;
|
||||
H{v}[i] = 0;
|
||||
}}
|
||||
"""
|
||||
|
||||
elif polarity > 0:
|
||||
bc_e = """
|
||||
if ({r} == s{r} - 1) {{
|
||||
E{r}[i] = -E{r}[i-2*di{r}];
|
||||
E{u}[i] = +E{u}[i-di{r}];
|
||||
E{v}[i] = +E{v}[i-di{r}];
|
||||
}} else if ({r} == s{r} - 2) {{
|
||||
E{r}[i] = 0;
|
||||
}}
|
||||
"""
|
||||
bc_h = """
|
||||
if ({r} == s{r} - 1) {{
|
||||
H{r}[i] = +H{r}[i-di{r}];
|
||||
H{u}[i] = -H{u}[i-2*di{r}];
|
||||
H{v}[i] = -H{v}[i-2*di{r}];
|
||||
}} else if ({r} == s{r} - 2) {{
|
||||
H{u}[i] = 0;
|
||||
H{v}[i] = 0;
|
||||
}}
|
||||
"""
|
||||
else:
|
||||
raise Exception()
|
||||
|
||||
replacements = {'r': r, 'u': uv[0], 'v': uv[1]}
|
||||
return [s.format(**replacements) for s in (bc_e, bc_h)]
|
||||
|
||||
|
||||
def cpml(direction: int,
|
||||
polarity: int,
|
||||
dt: float,
|
||||
thickness: int=8,
|
||||
epsilon_eff: float=1,
|
||||
) -> Dict:
|
||||
"""
|
||||
Generate source code for complex phase matched layer (cpml) absorbing boundaries.
|
||||
These are not full boundary conditions and require a conducting boundary to be added
|
||||
in the same direction as the pml.
|
||||
|
||||
:param direction: integer in range(3), corresponding to x, y, z directions.
|
||||
:param polarity: -1 or 1, corresponding to eg. -x or +x direction.
|
||||
:param dt: timestep used by the simulation
|
||||
:param thickness: Number of cells used by the pml (the grid is NOT expanded to add these cells). Default 8.
|
||||
:param epsilon_eff: Effective epsilon_r of the pml layer. Default 1.
|
||||
:return: Dict with entries 'E', 'H' (update equations for E and H) and 'psi_E', 'psi_H' (lists of str,
|
||||
specifying the field names of the cpml fields used in the E and H update steps. Eg.,
|
||||
Psi_xn_Ex for the complex Ex component for the x- pml.)
|
||||
"""
|
||||
if direction not in range(3):
|
||||
raise Exception('Invalid direction: {}'.format(direction))
|
||||
|
||||
if polarity not in (-1, 1):
|
||||
raise Exception('Invalid polarity: {}'.format(polarity))
|
||||
|
||||
if thickness <= 2:
|
||||
raise Exception('It would be wise to have a pml with 4+ cells of thickness')
|
||||
|
||||
if epsilon_eff <= 0:
|
||||
raise Exception('epsilon_eff must be positive')
|
||||
|
||||
m = (3.5, 1)
|
||||
sigma_max = 0.8 * (m[0] + 1) / numpy.sqrt(epsilon_eff)
|
||||
alpha_max = 0 # TODO: Decide what to do about non-zero alpha
|
||||
transverse = numpy.delete(range(3), direction)
|
||||
|
||||
r = 'xyz'[direction]
|
||||
np = 'nVp'[numpy.sign(polarity)+1]
|
||||
uv = ['xyz'[i] for i in transverse]
|
||||
|
||||
xe = numpy.arange(1, thickness+1, dtype=float)[::-1]
|
||||
xh = numpy.arange(1, thickness+1, dtype=float)[::-1]
|
||||
if polarity > 0:
|
||||
xe -= 0.5
|
||||
elif polarity < 0:
|
||||
xh -= 0.5
|
||||
|
||||
def par(x):
|
||||
sigma = ((x / thickness) ** m[0]) * sigma_max
|
||||
alpha = ((1 - x / thickness) ** m[1]) * alpha_max
|
||||
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]}
|
||||
|
||||
if polarity < 0:
|
||||
bounds_if = """
|
||||
if ( 0 < {r} && {r} < {th} + 1 ) {{
|
||||
const int ir = {r} - 1; // index into pml parameters
|
||||
const int ip = {v} + {u} * s{v} + ir * s{v} * s{u}; // linear index into Psi
|
||||
"""
|
||||
elif polarity > 0:
|
||||
bounds_if = """
|
||||
if ( (s{r} - 1) > {r} && {r} > (s{r} - 1) - ({th} + 1) ) {{
|
||||
const int ir = (s{r} - 1) - ({r} + 1); // index into pml parameters
|
||||
const int ip = {v} + {u} * s{v} + ir * s{v} * s{u}; // linear index into Psi
|
||||
"""
|
||||
else:
|
||||
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}]);
|
||||
|
||||
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]);
|
||||
|
||||
H{u}[i] {sh}= dt * Psi_{r}{np}_H{u}[ip];
|
||||
H{v}[i] {se}= dt * Psi_{r}{np}_H{v}[ip];
|
||||
}}
|
||||
"""
|
||||
|
||||
pml_data = {
|
||||
'E': (bounds_if + code_e).format(**vals),
|
||||
'H': (bounds_if + code_h).format(**vals),
|
||||
'psi_E': ['Psi_{r}{np}_E{u}'.format(**vals),
|
||||
'Psi_{r}{np}_E{v}'.format(**vals)],
|
||||
'psi_H': ['Psi_{r}{np}_H{u}'.format(**vals),
|
||||
'Psi_{r}{np}_H{v}'.format(**vals)],
|
||||
}
|
||||
|
||||
return pml_data
|
||||
@ -0,0 +1,84 @@
|
||||
{#
|
||||
/* Common code for E, H updates
|
||||
*
|
||||
* Template parameters:
|
||||
* ftype type name (e.g. float or double)
|
||||
* shape list of 3 ints specifying shape of fields
|
||||
*/
|
||||
#}
|
||||
|
||||
typedef {{ftype}} ftype;
|
||||
|
||||
|
||||
/*
|
||||
* Field size info
|
||||
*/
|
||||
const size_t sx = {{shape[0]}};
|
||||
const size_t sy = {{shape[1]}};
|
||||
const size_t sz = {{shape[2]}};
|
||||
const size_t field_size = sx * sy * sz;
|
||||
|
||||
//Since we use i to index into Ex[], Ey[], ... rather than E[], do nothing if
|
||||
// i is outside the bounds of Ex[].
|
||||
if (i >= field_size) {
|
||||
PYOPENCL_ELWISE_CONTINUE;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Array indexing
|
||||
*/
|
||||
// Given a linear index i and shape (sx, sy, sz), defines x, y, and z
|
||||
// as the 3D indices of the current element (i).
|
||||
// (ie, converts linear index [i] to field indices (x, y, z)
|
||||
const size_t x = i / (sz * sy);
|
||||
const size_t y = (i - x * sz * sy) / sz;
|
||||
const size_t z = (i - y * sz - x * sz * sy);
|
||||
|
||||
// Calculate linear index offsets corresponding to offsets in 3D
|
||||
// (ie, if E[i] <-> E(x, y, z), then E[i + diy] <-> E(x, y + 1, z)
|
||||
const size_t dix = sz * sy;
|
||||
const size_t diy = sz;
|
||||
const size_t diz = 1;
|
||||
|
||||
|
||||
/*
|
||||
* Pointer math
|
||||
*/
|
||||
//Pointer offsets into the components of a linearized vector-field
|
||||
// (eg. Hx = H + XX, where H and Hx are pointers)
|
||||
const size_t XX = 0;
|
||||
const size_t YY = field_size;
|
||||
const size_t ZZ = field_size * 2;
|
||||
|
||||
//Define pointers to vector components of each field (eg. Hx = H + XX)
|
||||
__global ftype *Ex = E + XX;
|
||||
__global ftype *Ey = E + YY;
|
||||
__global ftype *Ez = E + ZZ;
|
||||
|
||||
__global ftype *Hx = H + XX;
|
||||
__global ftype *Hy = H + YY;
|
||||
__global ftype *Hz = H + ZZ;
|
||||
|
||||
|
||||
/*
|
||||
* Implement periodic boundary conditions
|
||||
*
|
||||
* mx ([m]inus [x]) gives the index offset of the adjacent cell in the minus-x direction.
|
||||
* In the event that we start at x == 0, we actually want to wrap around and grab the cell
|
||||
* x_{-1} == (sx - 1) instead, ie. mx = (sx - 1) * dix .
|
||||
*
|
||||
* px ([p]lus [x]) gives the index offset of the adjacent cell in the plus-x direction.
|
||||
* In the event that we start at x == (sx - 1), we actually want to wrap around and grab
|
||||
* the cell x_{+1} == 0 instead, ie. px = -(sx - 1) * dix .
|
||||
*/
|
||||
{% for r in 'xyz' %}
|
||||
int m{{r}} = -di{{r}};
|
||||
int p{{r}} = +di{{r}};
|
||||
int wrap_{{r}} = (s{{r}} - 1) * di{{r}};
|
||||
if ( {{r}} == 0 ) {
|
||||
m{{r}} = wrap_{{r}};
|
||||
} else if ( {{r}} == s{{r}} - 1 ) {
|
||||
p{{r}} = -wrap_{{r}};
|
||||
}
|
||||
{% endfor %}
|
||||
@ -0,0 +1,78 @@
|
||||
/*
|
||||
* Update E-field, including any PMLs.
|
||||
*
|
||||
* Template parameters:
|
||||
* common_header: Rendered contents of common.cl
|
||||
* pmls: [('x', 'n'), ('z', 'p'),...] list of pml axes and polarities
|
||||
* pml_thickness: Number of cells (integer)
|
||||
*
|
||||
* OpenCL args:
|
||||
* E, H, dt, eps, [p{01}e{np}, Psi_{xyz}{np}_E]
|
||||
*/
|
||||
|
||||
{{common_header}}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
__global ftype *epsx = eps + XX;
|
||||
__global ftype *epsy = eps + YY;
|
||||
__global ftype *epsz = eps + ZZ;
|
||||
|
||||
{% if pmls -%}
|
||||
const int pml_thickness = {{pml_thickness}};
|
||||
{%- endif %}
|
||||
|
||||
/*
|
||||
* Precalclate derivatives
|
||||
*/
|
||||
ftype dHxy = Hx[i] - Hx[i + my];
|
||||
ftype dHxz = Hx[i] - Hx[i + mz];
|
||||
|
||||
ftype dHyx = Hy[i] - Hy[i + mx];
|
||||
ftype dHyz = Hy[i] - Hy[i + mz];
|
||||
|
||||
ftype dHzx = Hz[i] - Hz[i + mx];
|
||||
ftype dHzy = Hz[i] - Hz[i + my];
|
||||
|
||||
/*
|
||||
* PML Update
|
||||
*/
|
||||
// PML effects on E
|
||||
ftype pExi = 0;
|
||||
ftype pEyi = 0;
|
||||
ftype pEzi = 0;
|
||||
|
||||
{% for r, p in pmls -%}
|
||||
{%- set u, v = ['x', 'y', 'z'] | reject('equalto', r) -%}
|
||||
{%- set psi = 'Psi_' ~ r ~ p ~ '_E' -%}
|
||||
{%- if r != 'y' -%}
|
||||
{%- set se, sh = '-', '+' -%}
|
||||
{%- else -%}
|
||||
{%- set se, sh = '+', '-' -%}
|
||||
{%- endif -%}
|
||||
|
||||
{%- if p == 'n' %}
|
||||
|
||||
if ( {{r}} < pml_thickness ) {
|
||||
const size_t ir = {{r}}; // index into pml parameters
|
||||
|
||||
{%- elif p == 'p' %}
|
||||
|
||||
if ( s{{r}} > {{r}} && {{r}} >= s{{r}} - pml_thickness ) {
|
||||
const size_t ir = (s{{r}} - 1) - {{r}}; // index into pml parameters
|
||||
|
||||
{%- endif %}
|
||||
const size_t ip = {{v}} + {{u}} * s{{v}} + ir * s{{v}} * s{{u}}; // linear index into Psi
|
||||
{{psi ~ u}}[ip] = p0e{{p}}[ir] * {{psi ~ u}}[ip] + p1e{{p}}[ir] * dH{{v ~ r}};
|
||||
{{psi ~ v}}[ip] = p0e{{p}}[ir] * {{psi ~ v}}[ip] + p1e{{p}}[ir] * dH{{u ~ r}};
|
||||
pE{{u}}i {{se}}= {{psi ~ u}}[ip];
|
||||
pE{{v}}i {{sh}}= {{psi ~ v}}[ip];
|
||||
}
|
||||
{%- endfor %}
|
||||
|
||||
/*
|
||||
* Update E
|
||||
*/
|
||||
Ex[i] += dt / epsx[i] * (dHzy - dHyz + pExi);
|
||||
Ey[i] += dt / epsy[i] * (dHxz - dHzx + pEyi);
|
||||
Ez[i] += dt / epsz[i] * (dHyx - dHxy + pEzi);
|
||||
@ -0,0 +1,125 @@
|
||||
/*
|
||||
* Update H-field, including any PMLs.
|
||||
* Also precalculate values for poynting vector if necessary.
|
||||
*
|
||||
* Template parameters:
|
||||
* common_header: Rendered contents of common.cl
|
||||
* pmls: [('x', 'n'), ('z', 'p'),...] list of pml axes and polarities
|
||||
* pml_thickness: Number of cells (integer)
|
||||
* do_poynting: Whether to precalculate poynting vector components (boolean)
|
||||
*
|
||||
* OpenCL args:
|
||||
* E, H, dt, [p{01}h{np}, Psi_{xyz}{np}_H], [oS]
|
||||
*/
|
||||
|
||||
{{common_header}}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
{% if pmls -%}
|
||||
const int pml_thickness = {{pml_thickness}};
|
||||
{%- endif %}
|
||||
|
||||
/*
|
||||
* Precalculate derivatives
|
||||
*/
|
||||
ftype dExy = Ex[i + py] - Ex[i];
|
||||
ftype dExz = Ex[i + pz] - Ex[i];
|
||||
|
||||
ftype dEyx = Ey[i + px] - Ey[i];
|
||||
ftype dEyz = Ey[i + pz] - Ey[i];
|
||||
|
||||
ftype dEzx = Ez[i + px] - Ez[i];
|
||||
ftype dEzy = Ez[i + py] - Ez[i];
|
||||
|
||||
{%- if do_poynting %}
|
||||
|
||||
|
||||
/*
|
||||
* Precalculate averaged E
|
||||
*/
|
||||
ftype aExy = Ex[i + py] + Ex[i];
|
||||
ftype aExz = Ex[i + pz] + Ex[i];
|
||||
|
||||
ftype aEyx = Ey[i + px] + Ey[i];
|
||||
ftype aEyz = Ey[i + pz] + Ey[i];
|
||||
|
||||
ftype aEzx = Ez[i + px] + Ez[i];
|
||||
ftype aEzy = Ez[i + py] + Ez[i];
|
||||
{%- endif %}
|
||||
|
||||
|
||||
/*
|
||||
* PML Update
|
||||
*/
|
||||
// PML contributions to H
|
||||
ftype pHxi = 0;
|
||||
ftype pHyi = 0;
|
||||
ftype pHzi = 0;
|
||||
|
||||
{%- for r, p in pmls -%}
|
||||
{%- set u, v = ['x', 'y', 'z'] | reject('equalto', r) -%}
|
||||
{%- set psi = 'Psi_' ~ r ~ p ~ '_H' -%}
|
||||
{%- if r != 'y' -%}
|
||||
{%- set se, sh = '-', '+' -%}
|
||||
{%- else -%}
|
||||
{%- set se, sh = '+', '-' -%}
|
||||
{%- endif -%}
|
||||
|
||||
{%- if p == 'n' %}
|
||||
|
||||
if ( {{r}} < pml_thickness ) {
|
||||
const size_t ir = {{r}}; // index into pml parameters
|
||||
|
||||
{%- elif p == 'p' %}
|
||||
|
||||
if ( s{{r}} > {{r}} && {{r}} >= s{{r}} - pml_thickness ) {
|
||||
const size_t ir = (s{{r}} - 1) - {{r}}; // index into pml parameters
|
||||
|
||||
{%- endif %}
|
||||
const size_t ip = {{v}} + {{u}} * s{{v}} + ir * s{{v}} * s{{u}}; // linear index into Psi
|
||||
{{psi ~ u}}[ip] = p0h{{p}}[ir] * {{psi ~ u}}[ip] + p1h{{p}}[ir] * dE{{v ~ r}};
|
||||
{{psi ~ v}}[ip] = p0h{{p}}[ir] * {{psi ~ v}}[ip] + p1h{{p}}[ir] * dE{{u ~ r}};
|
||||
pH{{u}}i {{sh}}= {{psi ~ u}}[ip];
|
||||
pH{{v}}i {{se}}= {{psi ~ v}}[ip];
|
||||
}
|
||||
{%- endfor %}
|
||||
|
||||
/*
|
||||
* Update H
|
||||
*/
|
||||
{% if do_poynting -%}
|
||||
// Save old H for averaging
|
||||
ftype Hx_old = Hx[i];
|
||||
ftype Hy_old = Hy[i];
|
||||
ftype Hz_old = Hz[i];
|
||||
{%- endif %}
|
||||
|
||||
// H update equations
|
||||
Hx[i] -= dt * (dEzy - dEyz + pHxi);
|
||||
Hy[i] -= dt * (dExz - dEzx + pHyi);
|
||||
Hz[i] -= dt * (dEyx - dExy + pHzi);
|
||||
|
||||
{% if do_poynting -%}
|
||||
// Average H across timesteps
|
||||
ftype aHxt = Hx[i] + Hx_old;
|
||||
ftype aHyt = Hy[i] + Hy_old;
|
||||
ftype aHzt = Hz[i] + Hz_old;
|
||||
|
||||
/*
|
||||
* Calculate unscaled S components at H locations
|
||||
*/
|
||||
__global ftype *oSxy = oS + 0 * field_size;
|
||||
__global ftype *oSyz = oS + 1 * field_size;
|
||||
__global ftype *oSzx = oS + 2 * field_size;
|
||||
__global ftype *oSxz = oS + 3 * field_size;
|
||||
__global ftype *oSyx = oS + 4 * field_size;
|
||||
__global ftype *oSzy = oS + 5 * field_size;
|
||||
|
||||
oSxy[i] = aEyx * aHzt;
|
||||
oSxz[i] = -aEzx * aHyt;
|
||||
oSyz[i] = aEzy * aHxt;
|
||||
oSyx[i] = -aExy * aHzt;
|
||||
oSzx[i] = aExz * aHyt;
|
||||
oSzy[i] = -aEyz * aHxt;
|
||||
{%- endif -%}
|
||||
@ -0,0 +1,36 @@
|
||||
/*
|
||||
* Update E-field, including any PMLs.
|
||||
*
|
||||
* Template parameters:
|
||||
* common_header: Rendered contents of common.cl
|
||||
* pmls: [('x', 'n'), ('z', 'p'),...] list of pml axes and polarities
|
||||
* pml_thickness: Number of cells (integer)
|
||||
*
|
||||
* OpenCL args:
|
||||
* E, H, dt, S, oS
|
||||
*/
|
||||
|
||||
{{common_header}}
|
||||
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
/*
|
||||
* Calculate S from oS (pre-calculated components)
|
||||
*/
|
||||
__global ftype *Sx = S + XX;
|
||||
__global ftype *Sy = S + YY;
|
||||
__global ftype *Sz = S + ZZ;
|
||||
|
||||
// Use unscaled S components from H locations
|
||||
__global ftype *oSxy = oS + 0 * field_size;
|
||||
__global ftype *oSyz = oS + 1 * field_size;
|
||||
__global ftype *oSzx = oS + 2 * field_size;
|
||||
__global ftype *oSxz = oS + 3 * field_size;
|
||||
__global ftype *oSyx = oS + 4 * field_size;
|
||||
__global ftype *oSzy = oS + 5 * field_size;
|
||||
|
||||
ftype s_factor = dt * 0.125;
|
||||
Sx[i] = (oSxy[i] + oSxz[i] + oSxy[i + my] + oSxz[i + mz]) * s_factor;
|
||||
Sy[i] = (oSyz[i] + oSyx[i] + oSyz[i + mz] + oSyx[i + mx]) * s_factor;
|
||||
Sz[i] = (oSzx[i] + oSzy[i] + oSzx[i + mx] + oSzy[i + my]) * s_factor;
|
||||
@ -1,6 +1,8 @@
|
||||
numpy
|
||||
h5py
|
||||
pyopencl
|
||||
jinja2
|
||||
-e git+https://mpxd.net/gogs/jan/float_raster.git@release#egg=float_raster
|
||||
-e git+https://mpxd.net/gogs/jan/gridlock.git@release#egg=gridlock
|
||||
-e git+https://mpxd.net/gogs/jan/masque.git@release#egg=masque
|
||||
-e git+https://mpxd.net/gogs/jan/fdfd_tools.git@master#egg=fdfd_tools
|
||||
|
||||
Loading…
Reference in New Issue