diff --git a/opencl_fdtd/kernels/update_e.cl b/opencl_fdtd/kernels/update_e.cl index f44527d..9127181 100644 --- a/opencl_fdtd/kernels/update_e.cl +++ b/opencl_fdtd/kernels/update_e.cl @@ -5,9 +5,12 @@ * common_header: Rendered contents of common.cl * pmls: [{'axis': 'x', 'polarity': 'n', 'thickness': 8}, ...] list of pml dicts containing * axes, polarities, and thicknesses. + * uniform_dx: If grid is uniform, uniform_dx should be the grid spacing. + * Otherwise, uniform_dx should be False and [inv_dh{xyz}] arrays must be supplied as + * OpenCL args. * * OpenCL args: - * E, H, dt, eps, [p{01}e{np}, Psi_{xyz}{np}_E] + * E, H, dt, eps, [p{012}e{np}, Psi_{xyz}{np}_E], [inv_dh{xyz}] */ {{common_header}} @@ -19,17 +22,28 @@ __global ftype *epsy = eps + YY; __global ftype *epsz = eps + ZZ; +{%- if uniform_dx %} +ftype inv_dx = 1.0 / {{uniform_dx}}; +ftype inv_dy = 1.0 / {{uniform_dx}}; +ftype inv_dz = 1.0 / {{uniform_dx}}; +{%- else %} +ftype inv_dx = inv_dhx[x]; +ftype inv_dy = inv_dhy[y]; +ftype inv_dz = inv_dhz[z]; +{%- endif %} + + /* * Precalculate derivatives */ -ftype dHxy = Hx[i] - Hx[i + my]; -ftype dHxz = Hx[i] - Hx[i + mz]; +ftype dHxy = (Hx[i] - Hx[i + my]) * inv_dy; +ftype dHxz = (Hx[i] - Hx[i + mz]) * inv_dz; -ftype dHyx = Hy[i] - Hy[i + mx]; -ftype dHyz = Hy[i] - Hy[i + mz]; +ftype dHyx = (Hy[i] - Hy[i + mx]) * inv_dx; +ftype dHyz = (Hy[i] - Hy[i + mz]) * inv_dz; -ftype dHzx = Hz[i] - Hz[i + mx]; -ftype dHzy = Hz[i] - Hz[i + my]; +ftype dHzx = (Hz[i] - Hz[i + mx]) * inv_dx; +ftype dHzy = (Hz[i] - Hz[i + my]) * inv_dy; {% for bloch in bloch_boundaries -%} {%- set r = bloch['axis'] -%} diff --git a/opencl_fdtd/kernels/update_h.cl b/opencl_fdtd/kernels/update_h.cl index f1c19cb..ac75157 100644 --- a/opencl_fdtd/kernels/update_h.cl +++ b/opencl_fdtd/kernels/update_h.cl @@ -7,9 +7,12 @@ * pmls: [{'axis': 'x', 'polarity': 'n', 'thickness': 8}, ...] list of pml dicts containing * axes, polarities, and thicknesses. * do_poynting: Whether to precalculate poynting vector components (boolean) + * uniform_dx: If grid is uniform, uniform_dx should be the grid spacing. + * Otherwise, uniform_dx should be False and [inv_de{xyz}] arrays must be supplied as + * OpenCL args. * * OpenCL args: - * E, H, dt, [p{xyz}{01}h{np}, Psi_{xyz}{np}_H], [oS] + * E, H, dt, [inv_de{xyz}], [p{xyz}{01}h{np}, Psi_{xyz}{np}_H], [oS] */ {{common_header}} @@ -19,14 +22,25 @@ /* * Precalculate derivatives */ -ftype dExy = Ex[i + py] - Ex[i]; -ftype dExz = Ex[i + pz] - Ex[i]; +{%- if uniform_dx %} +ftype inv_dx = 1.0 / {{uniform_dx}}; +ftype inv_dy = 1.0 / {{uniform_dx}}; +ftype inv_dz = 1.0 / {{uniform_dx}}; +{%- else %} +ftype inv_dx = inv_dex[x]; +ftype inv_dy = inv_dey[y]; +ftype inv_dz = inv_dez[z]; +{%- endif %} -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]; +ftype dExy = (Ex[i + py] - Ex[i]) * inv_dy; +ftype dExz = (Ex[i + pz] - Ex[i]) * inv_dz; + +ftype dEyx = (Ey[i + px] - Ey[i]) * inv_dx; +ftype dEyz = (Ey[i + pz] - Ey[i]) * inv_dz; + +ftype dEzx = (Ez[i + px] - Ez[i]) * inv_dx; +ftype dEzy = (Ez[i + py] - Ez[i]) * inv_dy; {% for bloch in bloch_boundaries -%} diff --git a/opencl_fdtd/simulation.py b/opencl_fdtd/simulation.py index 0eaef5d..940725a 100644 --- a/opencl_fdtd/simulation.py +++ b/opencl_fdtd/simulation.py @@ -61,6 +61,7 @@ class Simulation(object): S = None # type: pyopencl.array.Array eps = None # type: pyopencl.array.Array dt = None # type: float + inv_dxes = None # type: List[pyopencl.array.Array] arg_type = None # type: numpy.float32 or numpy.float64 @@ -77,7 +78,8 @@ class Simulation(object): epsilon: List[numpy.ndarray], pmls: List[Dict[str, int or float]], bloch_boundaries: List[Dict[str, int or float]] = (), - dt: float = .99/numpy.sqrt(3), + dxes: List[List[numpy.ndarray]] or float = None, + dt: float = None, initial_fields: Dict[str, List[numpy.ndarray]] = None, context: pyopencl.Context = None, queue: pyopencl.CommandQueue = None, @@ -98,7 +100,7 @@ class Simulation(object): 'ln_R_per_layer': Desired (ln(R) / thickness) value. Default -1.6. 'm': Polynomial grading exponent. Default 3.5. 'ma': Exponent for alpha. Default 1. - :param dt: Time step. Default is .99/sqrt(3). + :param dt: Time step. Default is min(dxes) * .99/sqrt(3). :param initial_E: Initial E-field (default is 0 everywhere). Same format as epsilon. :param initial_H: Initial H-field (default is 0 everywhere). Same format as epsilon. :param context: pyOpenCL context. If not given, pyopencl.create_some_context(False) is called. @@ -124,7 +126,22 @@ class Simulation(object): self._create_context(context, queue) self._create_eps(epsilon) - if dt > .99/numpy.sqrt(3): + if dxes is None: + dxes = 1.0 + + if isinstance(dxes, (float, int)): + uniform_dx = dxes + min_dx = dxes + else: + uniform_dx = False + self.inv_dxes = [self._create_field(1 / dxn) for dxn in dxes[0] + dxes[1]] + min_dx = min(min(dxn) for dxn in dxes[0] + dxes[1]) + + max_dt = min_dx * .99 / numpy.sqrt(3) + + if dt is None: + self.dt = max_dt + elif dt > max_dt: warnings.warn('Warning: unstable dt: {}'.format(dt)) elif dt <= 0: raise Exception('Invalid dt: {}'.format(dt)) @@ -154,6 +171,10 @@ class Simulation(object): base_fields[ptr('E')] = self.E base_fields[ptr('H')] = self.H base_fields[ctype + ' dt'] = self.dt + if uniform_dx == 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): + base_fields[ptr(name)] = field eps_field = OrderedDict() eps_field[ptr('eps')] = self.eps @@ -178,6 +199,7 @@ class Simulation(object): 'pmls': pmls, 'do_poynting': do_poynting, 'bloch': bloch_boundaries, + 'uniform_dx': uniform_dx, } E_source = jinja_env.get_template('update_e.cl').render(**jinja_args) H_source = jinja_env.get_template('update_h.cl').render(**jinja_args)