Compare commits

..

No commits in common. "5b16295769d1b27f4cb267f511d561bee303a4b7" and "b703f1ee20576dcc899df67824586662f0e4a29e" have entirely different histories.

4 changed files with 72 additions and 131 deletions

View File

@ -1,7 +1,4 @@
from .simulation import (
Simulation as Simulation,
type_to_C as type_to_C,
)
from .simulation import Simulation, type_to_C
__author__ = 'Jan Petykiewicz'
__version__ = '0.4'

View File

@ -2,13 +2,12 @@
Class for constructing and holding the basic FDTD operations and fields
"""
from collections.abc import Callable, Sequence
from typing import Callable, Type, Sequence
from collections import OrderedDict
import numpy
from numpy.typing import NDArray
from numpy import floating, complexfloating
import jinja2
import warnings
import logging
import pyopencl
import pyopencl.array
@ -17,18 +16,13 @@ from pyopencl.elementwise import ElementwiseKernel
from meanas.fdmath import vec
logger = logging.getLogger(__name__)
__author__ = 'Jan Petykiewicz'
# Create jinja2 env on module load
jinja_env = jinja2.Environment(loader=jinja2.PackageLoader(__name__.split('.')[0], 'kernels'))
class FDTDError(Exception):
""" Custom exception for opencl_fdtd """
pass
class Simulation:
r"""
Constructs and holds the basic FDTD operations and related fields
@ -70,7 +64,7 @@ class Simulation:
dt: float
inv_dxes: list[pyopencl.array.Array]
arg_type: type
arg_type: Type
context: pyopencl.Context
queue: pyopencl.CommandQueue
@ -91,7 +85,7 @@ class Simulation:
initial_fields: dict[str, NDArray] | None = None,
context: pyopencl.Context | None = None,
queue: pyopencl.CommandQueue | None = None,
float_type: type = numpy.float32,
float_type: Type = numpy.float32,
do_poynting: bool = True,
do_fieldsrc: bool = False,
) -> None:
@ -140,7 +134,7 @@ class Simulation:
if dxes is None:
dxes = 1.0
if isinstance(dxes, float | int):
if isinstance(dxes, (float, int)):
uniform_dx = dxes
min_dx = dxes
else:
@ -149,14 +143,13 @@ class Simulation:
min_dx = min(min(dxn) for dxn in dxes[0] + dxes[1])
max_dt = min_dx * .99 / numpy.sqrt(3)
logger.info(f'{min_dx=}, {max_dt=}, {dt=}')
if dt is None:
self.dt = max_dt
elif dt > max_dt:
warnings.warn(f'Warning: unstable dt: {dt}', stacklevel=2)
warnings.warn(f'Warning: unstable dt: {dt}')
elif dt <= 0:
raise FDTDError(f'Invalid dt: {dt}')
raise Exception(f'Invalid dt: {dt}')
else:
self.dt = dt
@ -180,31 +173,28 @@ class Simulation:
def ptr(arg: str) -> str:
return ctype + ' *' + arg
base_fields = {
ptr('E'): self.E,
ptr('H'): self.H,
ctype + ' dt': self.dt,
}
base_fields = OrderedDict()
base_fields[ptr('E')] = self.E
base_fields[ptr('H')] = self.H
base_fields[ctype + ' dt'] = self.dt
if uniform_dx is False:
inv_dx_names = ['inv_d' + eh + r for eh in 'eh' for r in 'xyz']
for name, field in zip(inv_dx_names, self.inv_dxes, strict=True):
for name, field in zip(inv_dx_names, self.inv_dxes):
base_fields[ptr(name)] = field
eps_field = {ptr('eps'): self.eps}
eps_field = OrderedDict()
eps_field[ptr('eps')] = self.eps
if bloch_boundaries:
base_fields |= {
ptr('F'): self.F,
ptr('G'): self.G,
}
base_fields[ptr('F')] = self.F
base_fields[ptr('G')] = self.G
bloch_fields = {
ptr('E'): self.F,
ptr('H'): self.G,
ctype + ' dt': self.dt,
ptr('F'): self.E,
ptr('G'): self.H,
}
bloch_fields = OrderedDict()
bloch_fields[ptr('E')] = self.F
bloch_fields[ptr('H')] = self.G
bloch_fields[ctype + ' dt'] = self.dt
bloch_fields[ptr('F')] = self.E
bloch_fields[ptr('G')] = self.H
common_source = jinja_env.get_template('common.cl').render(
ftype=ctype,
@ -226,8 +216,8 @@ class Simulation:
if bloch_boundaries:
bloch_args = jinja_args.copy()
bloch_args['do_poynting'] = False
bloch_args['bloch'] = [{
'axis': b['axis'],
bloch_args['bloch'] = [
{'axis': b['axis'],
'real': b['imag'],
'imag': b['real'],
}
@ -237,7 +227,7 @@ class Simulation:
self.sources['F'] = F_source
self.sources['G'] = G_source
S_fields = {}
S_fields = OrderedDict()
if do_poynting:
self.S = pyopencl.array.zeros_like(self.E)
S_fields[ptr('S')] = self.S
@ -247,7 +237,7 @@ class Simulation:
S_fields[ptr('S0')] = self.S0
S_fields[ptr('S1')] = self.S1
J_fields = {}
J_fields = OrderedDict()
if do_fieldsrc:
J_source = jinja_env.get_template('update_j.cl').render(**jinja_args)
self.sources['J'] = J_source
@ -257,36 +247,37 @@ class Simulation:
J_fields[ptr('Jr')] = self.Jr
J_fields[ptr('Ji')] = self.Ji
#
# PML
#
'''
PML
'''
pml_e_fields, pml_h_fields = self._create_pmls(pmls)
if bloch_boundaries:
pml_f_fields, pml_g_fields = self._create_pmls(pmls)
#
# Create operations
#
'''
Create operations
'''
self.update_E = self._create_operation(E_source, (base_fields, eps_field, pml_e_fields))
self.update_H = self._create_operation(H_source, (base_fields, pml_h_fields, S_fields))
if bloch_boundaries:
self.update_F = self._create_operation(F_source, (bloch_fields, eps_field, pml_f_fields))
self.update_G = self._create_operation(G_source, (bloch_fields, pml_g_fields))
if do_fieldsrc:
args = base_fields | J_fields
args = OrderedDict()
[args.update(d) for d in (base_fields, J_fields)]
var_args = [ctype + ' ' + v for v in 'cs'] + ['uint ' + r + m for r in 'xyz' for m in ('min', 'max')]
update = ElementwiseKernel(self.context, operation=J_source,
arguments=', '.join(list(args.keys()) + var_args))
self.update_J = lambda e, *a: update(*args.values(), *a, wait_for=e)
def _create_pmls(self, pmls: Sequence[dict[str, float]]) -> tuple[dict[str, pyopencl.array.Array], ...]:
def _create_pmls(self, pmls):
ctype = type_to_C(self.arg_type)
def ptr(arg: str) -> str:
return ctype + ' *' + arg
pml_e_fields = {}
pml_h_fields = {}
pml_e_fields = OrderedDict()
pml_h_fields = OrderedDict()
for pml in pmls:
a = 'xyz'.find(pml['axis'])
@ -294,9 +285,7 @@ class Simulation:
kappa_max = numpy.sqrt(pml['mu_eff'] * pml['epsilon_eff'])
alpha_max = pml['cfs_alpha']
print(sigma_max, kappa_max, alpha_max, pml['thickness'], self.dt)
def par(x, pml=pml, sigma_max=sigma_max, kappa_max=kappa_max, alpha_max=alpha_max): # noqa: ANN001, ANN202
def par(x):
scaling = (x / pml['thickness']) ** pml['m']
sigma = scaling * sigma_max
kappa = 1 + scaling * (kappa_max - 1)
@ -312,13 +301,8 @@ class Simulation:
elif pml['polarity'] == 'n':
xh -= 0.5
logger.debug(f'{pml=}')
logger.debug(f'{xe=}')
logger.debug(f'{xh=}')
logger.debug(f'{par(xe)=}')
logger.debug(f'{par(xh)=}')
pml_p_names = [['p' + pml['axis'] + i + eh + pml['polarity'] for i in '012'] for eh in 'eh']
for name_e, name_h, pe, ph in zip(pml_p_names[0], pml_p_names[1], par(xe), par(xh), strict=True):
for name_e, name_h, pe, ph in zip(pml_p_names[0], pml_p_names[1], par(xe), par(xh)):
pml_e_fields[ptr(name_e)] = pyopencl.array.to_device(self.queue, pe)
pml_h_fields[ptr(name_h)] = pyopencl.array.to_device(self.queue, ph)
@ -329,13 +313,13 @@ class Simulation:
psi_shape = list(self.shape)
psi_shape[a] = pml['thickness']
for ne, nh in zip(*psi_names, strict=True):
for ne, nh in zip(*psi_names):
pml_e_fields[ptr(ne)] = pyopencl.array.zeros(self.queue, tuple(psi_shape), dtype=self.arg_type)
pml_h_fields[ptr(nh)] = pyopencl.array.zeros(self.queue, tuple(psi_shape), dtype=self.arg_type)
return pml_e_fields, pml_h_fields
def _create_operation(self, source: str, args_fields: Sequence[dict[str, pyopencl.array.Array]]) -> Callable[..., pyopencl.Event]:
args = {}
def _create_operation(self, source, args_fields) -> Callable[..., pyopencl.Event]:
args = OrderedDict()
for d in args_fields:
args.update(d)
update = ElementwiseKernel(
@ -350,30 +334,38 @@ class Simulation:
context: pyopencl.Context | None = None,
queue: pyopencl.CommandQueue | None = None,
) -> None:
self.context = context or pyopencl.create_some_context()
self.queue = queue or pyopencl.CommandQueue(self.context)
if context is None:
self.context = pyopencl.create_some_context()
else:
self.context = context
if queue is None:
self.queue = pyopencl.CommandQueue(self.context)
else:
self.queue = queue
def _create_eps(self, epsilon: NDArray) -> pyopencl.array.Array:
if len(epsilon) != 3:
raise FDTDError('Epsilon must be a list with length of 3')
if not all(e.shape == epsilon[0].shape for e in epsilon[1:]):
raise FDTDError('All epsilon grids must have the same shape. Shapes are {}', [e.shape for e in epsilon])
raise Exception('Epsilon must be a list with length of 3')
if not all((e.shape == epsilon[0].shape for e in epsilon[1:])):
raise Exception('All epsilon grids must have the same shape. Shapes are {}', [e.shape for e in epsilon])
if not epsilon[0].shape == self.shape:
raise FDTDError(f'Epsilon shape mismatch. Expected {self.shape}, got {epsilon[0].shape}')
raise Exception(f'Epsilon shape mismatch. Expected {self.shape}, got {epsilon[0].shape}')
self.eps = pyopencl.array.to_device(self.queue, vec(epsilon).astype(self.arg_type))
def _create_field(self, initial_value: NDArray | None = None) -> pyopencl.array.Array:
if initial_value is None:
return pyopencl.array.zeros_like(self.eps)
else:
if len(initial_value) != 3:
raise FDTDError('Initial field value must be a list of length 3')
if not all(f.shape == self.shape for f in initial_value):
raise FDTDError('Initial field list elements must have same shape as epsilon elements')
Exception('Initial field value must be a list of length 3')
if not all((f.shape == self.shape for f in initial_value)):
Exception('Initial field list elements must have same shape as epsilon elements')
return pyopencl.array.to_device(self.queue, vec(initial_value).astype(self.arg_type))
def type_to_C(
float_type: type,
float_type: Type,
) -> str:
"""
Returns a string corresponding to the C equivalent of a numpy type.
@ -393,7 +385,7 @@ def type_to_C(
numpy.complex128: 'cdouble_t',
}
if float_type not in types:
raise FDTDError(f'Unsupported type: {float_type}')
raise Exception(f'Unsupported type: {float_type}')
return types[float_type]

View File

@ -323,9 +323,8 @@ class JdotE_3DUniformDX(unittest.TestCase):
self.epsilon = numpy.full(shape, 4, dtype=float)
self.epsilon[self.src_mask] = 2
rng = numpy.random.default_rng()
e = rng.integers(-128, 128 + 1, size=shape, dtype=float)
h = rng.integers(-128, 128 + 1, size=shape, dtype=float)
e = numpy.random.randint(-128, 128 + 1, size=shape).astype(float)
h = numpy.random.randint(-128, 128 + 1, size=shape).astype(float)
self.es = [e]
self.hs = [h]
self.ss = []

View File

@ -32,10 +32,10 @@ classifiers = [
"License :: OSI Approved :: GNU Affero General Public License v3",
"Topic :: Scientific/Engineering",
]
requires-python = ">=3.11"
requires-python = ">=3.8"
dynamic = ["version"]
dependencies = [
"numpy>=1.26",
"numpy~=1.21",
"pyopencl",
"jinja2",
"meanas>=0.3",
@ -44,50 +44,3 @@ dependencies = [
[tool.hatch.version]
path = "opencl_fdtd/__init__.py"
[tool.ruff]
exclude = [
".git",
"dist",
]
line-length = 145
indent-width = 4
lint.dummy-variable-rgx = "^(_+|(_+[a-zA-Z0-9_]*[a-zA-Z0-9]+?))$"
lint.select = [
"NPY", "E", "F", "W", "B", "ANN", "UP", "SLOT", "SIM", "LOG",
"C4", "ISC", "PIE", "PT", "RET", "TCH", "PTH", "INT",
"ARG", "PL", "R", "TRY",
"G010", "G101", "G201", "G202",
"Q002", "Q003", "Q004",
]
lint.ignore = [
#"ANN001", # No annotation
"ANN002", # *args
"ANN003", # **kwargs
"ANN401", # Any
"ANN101", # self: Self
"SIM108", # single-line if / else assignment
"RET504", # x=y+z; return x
"PIE790", # unnecessary pass
"ISC003", # non-implicit string concatenation
"C408", # dict(x=y) instead of {'x': y}
"PLR09", # Too many xxx
"PLR2004", # magic number
"PLC0414", # import x as x
"TRY003", # Long exception message
]
[[tool.mypy.overrides]]
module = [
"scipy",
"scipy.optimize",
"scipy.linalg",
"scipy.sparse",
"scipy.sparse.linalg",
"pyopencl",
"pyopencl.array",
"pyopencl.elementwise",
"pyopencl.reduction",
]
ignore_missing_imports = true