From 4af25d5de1dcfa39bfcdeb191f4274db392a486c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 13 Dec 2016 13:50:11 -0600 Subject: [PATCH 1/8] Reorganize photonics code, start on SKIE --- examples/find-photonic-mode-sk.py | 123 +++ .../pde/{maxwell.py => maxwell/__init__.py} | 14 +- pytential/symbolic/pde/maxwell/fiber.py | 857 ++++++++++++++++++ pytential/symbolic/pde/scalar.py | 626 ------------- 4 files changed, 981 insertions(+), 639 deletions(-) create mode 100644 examples/find-photonic-mode-sk.py rename pytential/symbolic/pde/{maxwell.py => maxwell/__init__.py} (97%) create mode 100644 pytential/symbolic/pde/maxwell/fiber.py diff --git a/examples/find-photonic-mode-sk.py b/examples/find-photonic-mode-sk.py new file mode 100644 index 00000000..c44a6ce8 --- /dev/null +++ b/examples/find-photonic-mode-sk.py @@ -0,0 +1,123 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = """ +Copyright (C) 2015 Shidong Jiang, Andreas Kloeckner +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np +import numpy.linalg as la +import pyopencl as cl +from pytential import sym +from functools import partial + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl + as pytest_generate_tests) + + +def find_mode(): + import warnings + warnings.simplefilter("error", np.ComplexWarning) + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + k0 = 1.4447 + k1 = k0*1.02 + beta_sym = sym.var("beta") + + from pytential.symbolic.pde.maxwell.fiber import \ + SecondKindInfZMuellerOperator + + pde_op = SecondKindInfZMuellerOperator( + k_vacuum=1, + interfaces=((0, 1, sym.DEFAULT_SOURCE),), + domain_k_exprs=(k0, k1), + beta=beta_sym, + use_l2_weighting=True) + + jm_sym = pde_op.make_unknown("") + op = pde_op.operator(jm_sym) + 1/0 + + # {{{ discretization setup + + from meshmode.mesh.generation import ellipse, make_curve_mesh + curve_f = partial(ellipse, 1) + + target_order = 7 + qbx_order = 4 + nelements = 30 + + from meshmode.mesh.processing import affine_map + mesh = make_curve_mesh(curve_f, + np.linspace(0, 1, nelements+1), + target_order) + lambda_ = 1.55 + circle_radius = 3.4*2*np.pi/lambda_ + mesh = affine_map(mesh, A=circle_radius*np.eye(2)) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + from pytential.qbx import QBXLayerPotentialSource + density_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + qbx = QBXLayerPotentialSource(density_discr, 4*target_order, + qbx_order, + # Don't use FMM for now + fmm_order=False) + + # }}} + + x_vec = np.random.randn(len(u_sym)*density_discr.nnodes) + y_vec = np.random.randn(len(u_sym)*density_discr.nnodes) + + def muller_solve_func(beta): + from pytential.symbolic.execution import build_matrix + mat = build_matrix( + queue, qbx, op, u_sym, + context={"beta": beta}).get() + + return 1/x_vec.dot(la.solve(mat, y_vec)) + + starting_guesses = (1+0j)*( + k0 + + (k1-k0) * np.random.rand(3)) + + from pytential.muller import muller + beta, niter = muller(muller_solve_func, z_start=starting_guesses) + print("beta") + + +if __name__ == "__main__": + import sys + if sys.argv[1:]: + + exec(sys.argv[1]) + else: + find_mode() + +# vim: fdm=marker diff --git a/pytential/symbolic/pde/maxwell.py b/pytential/symbolic/pde/maxwell/__init__.py similarity index 97% rename from pytential/symbolic/pde/maxwell.py rename to pytential/symbolic/pde/maxwell/__init__.py index ea96d3ef..68cceb5e 100644 --- a/pytential/symbolic/pde/maxwell.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -25,19 +25,7 @@ THE SOFTWARE. from six.moves import range, zip import numpy as np # noqa from pymbolic.primitives import Variable -from pytential.symbolic.primitives import ( - cse, - xyz_to_tangential, tangential_to_xyz, n_cross, n_dot, - make_normal, make_tangent, make_vector_field, project_to_tangential, - - IterativeInverse, - - S, Sp, nxcurl_S, grad_S, curl_S_volume, surf_grad_S, - IntGdSource, - surface_laplacian_S_squared, S_surface_laplacian_S, - - Integral, LineIntegral, Mean, Ones, - ) +from pytential import sym # {{{ MFIE diff --git a/pytential/symbolic/pde/maxwell/fiber.py b/pytential/symbolic/pde/maxwell/fiber.py new file mode 100644 index 00000000..cd11ca63 --- /dev/null +++ b/pytential/symbolic/pde/maxwell/fiber.py @@ -0,0 +1,857 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2013-2016 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + +__doc__ = """ +Second-Kind Fiberoptic +^^^^^^^^^^^^^^^^^^^^^^ + +.. autoclass:: SecondKindInfZMuellerOperator + +2D Dielectric (old-style) +^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. autoclass:: DielectricSRep2DBoundaryOperator +.. autoclass:: DielectricSDRep2DBoundaryOperator +""" + +import numpy as np +from collections import namedtuple +from six.moves import range + +from pytential import sym +from pytential.symbolic.primitives import ( + cse, + sqrt_jac_q_weight, QWeight, area_element) +from pytools import memoize_method +from pytential.symbolic.pde.scalar import L2WeightedPDEOperator + + +# {{{ second-kind infinite-z Mueller operator + +class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): + """ + Second-kind IE representation by + `Lai and Jiang `_. + """ + + def __init__(self, k_vacuum, domain_k_exprs, beta, + interfaces, use_l2_weighting=None): + """ + :attr k_vacuum: A symbolic expression for the wave number in vacuum. + May be a string, which will be interpreted as a variable name. + :attr interfaces: a tuple of tuples + ``(outer_domain, inner_domain, interface_id)``, + where *outer_domain* and *inner_domain* are indices into + *domain_k_names*, + and *interface_id* is a symbolic name for the discretization of the + interface. 'outer' designates the side of the interface to which + the normal points. + :attr domain_k_exprs: a tuple of variable names of the Helmholtz + parameter *k*, to be used inside each part of the source geometry. + May also be a tuple of strings, which will be transformed into + variable references of the corresponding names. + :attr beta: A symbolic expression for the wave number in the :math:`z` + direction. May be a string, which will be interpreted as a variable + name. + """ + + self.interfaces = interfaces + + if isinstance(beta, str): + beta = sym.var(beta) + self.beta = sym.cse(beta, "beta") + + if isinstance(k_vacuum, str): + k_vacuum = sym.var(k_vacuum) + self.k_vacuum = sym.cse(k_vacuum, "k_vac") + + self.domain_k_exprs = [ + sym.var(k_expr) + if isinstance(k_expr, str) + else sym.cse(k_expr, "k%d" % idom) + for idom, k_expr in enumerate(domain_k_exprs)] + del domain_k_exprs + + self.domain_K_exprs = [ + sym.cse((k_expr**2-beta**2)**0.5, "K%d" % i) + for i, k_expr in enumerate(self.domain_k_exprs)] + + from sumpy.kernel import HelmholtzKernel + self.kernel = HelmholtzKernel(2, allow_evanescent=True) + + def make_unknown(self, prefix): + from pytools import reverse_dictionary + index_to_name_and_intf = reverse_dictionary(self.unknown_index) + + from pytools.obj_array import make_obj_array + return make_obj_array([ + sym.var(prefix + + index_to_name_and_intf[i][0] + + str(index_to_name_and_intf[i][1])) + for i in range(len(index_to_name_and_intf)) + ]) + + @property + @memoize_method + def unknown_index(self): + result = {} + for i in range(len(self.interfaces)): + for current in ["j", "m"]: + for cur_dir in ["z", "t"]: + result[(current + cur_dir, i)] = len(result) + + return result + + def representation(self, unknown): + raise NotImplementedError() + + def operator(self, unknown): + result = np.zeros(4*len(self.interfaces)) + + unk_idx = self.unknown_index + + for i in range(len(self.interfaces)): + idx_jt = unk_idx["jt", i] + idx_jz = unk_idx["jz", i] + idx_mt = unk_idx["mt", i] + idx_mz = unk_idx["mz", i] + + phi1 = jt = unknown[idx_jt] + phi2 = jz = unknown[idx_jz] + phi3 = mt = unknown[idx_mt] + phi4 = mz = unknown[idx_mz] + + ne = self.beta/self.k_vacuum + + dom0_idx, dom1_idx, where = self.interfaces[i] + dom_indices = [dom0_idx, dom1_idx] + + def S(dom, density): # noqa + return sym.S( + self.kernel, density, + k=self.domain_K_exprs[dom_indices[dom]]) + + def D(dom, density): # noqa + return sym.D( + self.kernel, density, + k=self.domain_K_exprs[dom_indices[dom]]) + + def T(dom, density): # noqa + 1/0 + + def Tt(dom, density): # noqa + 1/0 + + def Sn(dom, density): # noqa + 1/0 + + def St(dom, density): # noqa + 1/0 + + tangent = (sym.pseudoscalar(2, 1, where) + / sym.area_element(2, 1, where)) + normal = sym.normal(2, 1, where) + print(sym.pretty( + )) + 1/0 + + n0 = 1 # FIXME + n1 = 2 # FIXME + + a11 = n0**2 * D(0, phi1) - n1**2 * D(1, phi1) + a22 = -n0**2 * Sn(0, phi2) + n1**2 * Sn(1, phi2) + a33 = D(0, phi3)-D(1, phi3) + a44 = -Sn(0, phi4) + Sn(1, phi4) + + a21 = -1j * ne * ( + n0**2 * tangent.scalar_product( + S(0, normal * phi1)) + - n1**2 * tangent.scalar_product( + S(1, normal * phi1))) + a43 = -1j * ne * ( + tangent.scalar_product( + S(0, normal * phi3)) + - tangent.scalar_product( + S(1, normal * phi3))) + a13 = ne*(T(0, phi3) - T(1, phi3)) + a24 = ne*(St(0, phi4) - St(1, phi4)) + a14 = 1j*( + (n0**2 - ne**2) * S(0, phi4) + - (n1**2 - ne**2) * S(1, phi4) + ) + + a23 = ( + 1j * (Tt(0, phi3) - Tt(1, phi3)) + - 1j * ne * ( + n0**2 * tangent.scalar_product( + S(0, normal * phi3)) + - n1**2 * tangent.scalar_product( + S(1, normal * phi3)))) + a31 = -a13 + a32 = -a14 + a41 = -a23 + a42 = -a24 + + d1 = (n0**2 + n1**2)*phi1 + d2 = (n0**2 + n1**2)*phi2 + d3 = phi3 + d4 = phi4 + + result[idx_jt] += d1 + a11 + 000 + a13 + a14 + result[idx_jz] += d2 + a21 + a22 + a23 + a24 + result[idx_mt] += d3 + a31 + a32 + a33 + 0 + result[idx_mz] += d4 + a41 + a42 + a43 + a44 + + # TODO: L2 weighting + +# }}} + + +# {{{ old-style fiber + +class Dielectric2DBoundaryOperatorBase(L2WeightedPDEOperator): + r""" + Solves the following system of BVPs on :math:`\mathbb{R}^2`, in which + a disjoint family of domains :math:`\Omega_i` is embedded: + + .. math:: + + \triangle E + (k_i^2-\beta^2) E = 0\quad \text{on $\Omega_i$}\\ + \triangle H + (k_i^2-\beta^2) H = 0\quad \text{on $\Omega_i$}\\ + [H] = 0 \text{ on $\partial \Omega_i$},\quad + [E] = 0 \text{ on $\partial \Omega_i$}\\ + \left[ \frac{k_0}{k^2-\beta^2} \partial_{\hat n}H\right] = 0 + \quad\text{on $\partial \Omega_i$},\quad\\ + \left[ \frac{k_0}{k^2-\beta^2} \partial_{\hat n}E\right] = 0 + \quad\text{on $\partial \Omega_i$} + + :math:`E` and :math:`H` are assumed to be of the form + + .. math:: + + E(x,y,z,t)=E(x,y)e^{i(\beta z-\omega t) + H(x,y,z,t)=H(x,y)e^{i(\beta z-\omega t) + + where :math:`[\cdot]` denotes the jump across an interface, and :math:`k` + (without an index) denotes the value of :math:`k` on either side of the + interface, for the purpose of computing the jump. :math:`\hat n` denotes + the unit normal of the interface. + + .. automethod:: make_unknown + .. automethod:: representation_outer + .. automethod:: representation_inner + .. automethod:: operator + """ + + field_kind_e = 0 + field_kind_h = 1 + field_kinds = [field_kind_e, field_kind_h] + + side_in = 0 + side_out = 1 + sides = [side_in, side_out] + side_to_sign = { + side_in: -1, + side_out: 1, + } + + dir_none = 0 + dir_normal = 1 + dir_tangential = 2 + + BCTermDescriptor = namedtuple("BCDescriptor", + "i_interface direction field_kind coeff_inner coeff_outer".split()) + + # {{{ constructor + + def __init__(self, mode, k_vacuum, domain_k_exprs, beta, + interfaces, use_l2_weighting=None): + """ + :attr mode: one of 'te', 'tm', 'tem' + :attr k_vacuum: A symbolic expression for the wave number in vacuum. + May be a string, which will be interpreted as a variable name. + :attr interfaces: a tuple of tuples + ``(outer_domain, inner_domain, interface_id)``, + where *outer_domain* and *inner_domain* are indices into + *domain_k_names*, + and *interface_id* is a symbolic name for the discretization of the + interface. 'outer' designates the side of the interface to which + the normal points. + :attr domain_k_exprs: a tuple of variable names of the Helmholtz + parameter *k*, to be used inside each part of the source geometry. + May also be a tuple of strings, which will be transformed into + variable references of the corresponding names. + :attr beta: A symbolic expression for the wave number in the :math:`z` + direction. May be a string, which will be interpreted as a variable + name. + """ + + if use_l2_weighting is None: + use_l2_weighting = False + + super(Dielectric2DBoundaryOperatorBase, self).__init__( + use_l2_weighting=use_l2_weighting) + + if mode == "te": + self.ez_enabled = False + self.hz_enabled = True + elif mode == "tm": + self.ez_enabled = True + self.hz_enabled = False + elif mode == "tem": + self.ez_enabled = True + self.hz_enabled = True + else: + raise ValueError("invalid mode '%s'" % mode) + + self.interfaces = interfaces + + fk_e = self.field_kind_e + fk_h = self.field_kind_h + + dir_none = self.dir_none + dir_normal = self.dir_normal + dir_tangential = self.dir_tangential + + if isinstance(beta, str): + beta = sym.var(beta) + beta = sym.cse(beta, "beta") + + if isinstance(k_vacuum, str): + k_vacuum = sym.var(k_vacuum) + k_vacuum = sym.cse(k_vacuum, "k_vac") + + self.domain_k_exprs = [ + sym.var(k_expr) + if isinstance(k_expr, str) + else sym.cse(k_expr, "k%d" % idom) + for idom, k_expr in enumerate(domain_k_exprs)] + del domain_k_exprs + + # Note the case of k/K! + # "K" is the 2D Helmholtz parameter. + # "k" is the 3D Helmholtz parameter. + + self.domain_K_exprs = [ + sym.cse((k_expr**2-beta**2)**0.5, "K%d" % i) + for i, k_expr in enumerate(self.domain_k_exprs)] + + from sumpy.kernel import HelmholtzKernel + self.kernel = HelmholtzKernel(2, allow_evanescent=True) + + # {{{ build bc list + + # list of tuples, where each tuple consists of BCTermDescriptor instances + + all_bcs = [] + for i_interface, (outer_domain, inner_domain, _) in ( + enumerate(self.interfaces)): + k_outer = self.domain_k_exprs[outer_domain] + k_inner = self.domain_k_exprs[inner_domain] + + all_bcs += [ + ( # [E] = 0 + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_none, + field_kind=fk_e, + coeff_outer=1, + coeff_inner=-1), + ), + ( # [H] = 0 + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_none, + field_kind=fk_h, + coeff_outer=1, + coeff_inner=-1), + ), + ( + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_tangential, + field_kind=fk_e, + coeff_outer=beta/(k_outer**2-beta**2), + coeff_inner=-beta/(k_inner**2-beta**2)), + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_normal, + field_kind=fk_h, + coeff_outer=sym.cse(-k_vacuum/(k_outer**2-beta**2)), + coeff_inner=sym.cse(k_vacuum/(k_inner**2-beta**2))), + ), + ( + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_tangential, + field_kind=fk_h, + coeff_outer=beta/(k_outer**2-beta**2), + coeff_inner=-beta/(k_inner**2-beta**2)), + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_normal, + field_kind=fk_e, + coeff_outer=sym.cse( + (k_outer**2/k_vacuum)/(k_outer**2-beta**2)), + coeff_inner=sym.cse( + -(k_inner**2/k_vacuum) + / (k_inner**2-beta**2))) + ), + ] + + del k_outer + del k_inner + + self.bcs = [] + for bc in all_bcs: + any_significant_e = any( + term.field_kind == fk_e + and term.direction in [dir_normal, dir_none] + for term in bc) + any_significant_h = any( + term.field_kind == fk_h + and term.direction in [dir_normal, dir_none] + for term in bc) + is_necessary = ( + (self.ez_enabled and any_significant_e) + or + (self.hz_enabled and any_significant_h)) + + # Only keep tangential modes for TEM. Otherwise, + # no jump in H already implies jump condition on + # tangential derivative. + is_tem = self.ez_enabled and self.hz_enabled + terms = tuple( + term + for term in bc + if term.direction != dir_tangential + or is_tem) + + if is_necessary: + self.bcs.append(terms) + + assert (len(all_bcs) + * (int(self.ez_enabled) + int(self.hz_enabled)) // 2 + == len(self.bcs)) + + # }}} + + # }}} + + def is_field_present(self, field_kind): + return ( + (field_kind == self.field_kind_e and self.ez_enabled) + or + (field_kind == self.field_kind_h and self.hz_enabled)) + + def make_unknown(self, name): + num_densities = ( + 2 + * (int(self.ez_enabled) + int(self.hz_enabled)) + * len(self.interfaces)) + + assert num_densities == len(self.bcs) + + return sym.make_sym_vector(name, num_densities) + + def bc_term_to_operator_contrib(self, term, side, raw_potential_op, + density, discrete): + potential_op = raw_potential_op + + side_sign = self.side_to_sign[side] + + domain_outer, domain_inner, interface_id = \ + self.interfaces[term.i_interface] + if side == self.side_in: + K_expr = self.domain_K_exprs[domain_inner] # noqa + bc_coeff = term.coeff_inner + elif side == self.side_out: + K_expr = self.domain_K_exprs[domain_outer] # noqa + bc_coeff = term.coeff_outer + else: + raise ValueError("invalid value of 'side'") + + potential_op = potential_op( + self.kernel, density, source=interface_id, + k=K_expr) + + if term.direction == self.dir_none: + if raw_potential_op is sym.S: + jump_term = 0 + elif raw_potential_op is sym.D: + jump_term = (side_sign*0.5) * discrete + else: + assert False, raw_potential_op + elif term.direction == self.dir_normal: + potential_op = sym.normal_derivative( + potential_op, interface_id) + + if raw_potential_op is sym.S: + # S' + jump_term = (-side_sign*0.5) * discrete + elif raw_potential_op is sym.D: + jump_term = 0 + else: + assert False, raw_potential_op + + elif term.direction == self.dir_tangential: + potential_op = sym.tangential_derivative( + raw_potential_op( + self.kernel, density, source=interface_id, + k=K_expr, qbx_forced_limit=side_sign), + interface_id).a.as_scalar() + + # Some of these may have jumps, but QBX does the dirty + # work here by directly computing the limit. + jump_term = 0 + + else: + raise ValueError("invalid direction") + + potential_op = ( + jump_term + + self.get_sqrt_weight(interface_id)*potential_op) + + del jump_term + + contrib = bc_coeff * potential_op + + if (raw_potential_op is sym.D + and term.direction == self.dir_normal): + # FIXME The hypersingular part should perhaps be + # treated specially to avoid cancellation. + pass + + return contrib + + +# {{{ single-layer representation + +class DielectricSRep2DBoundaryOperator(Dielectric2DBoundaryOperatorBase): + def _structured_unknown(self, unknown, with_l2_weights): + """ + :arg with_l2_weights: If True, return the 'bare' unknowns + that do not have the :math:`L^2` weights divided out. + Note: Those unknowns should *not* be interpreted as + point values of a density. + :returns: an array of unknowns, with the following index axes: + ``[side, field_kind, i_interface]``, where + ``side`` is o for the outside part and i for the interior part, + ``field_kind`` is 0 for the E-field and 1 for the H-field part, + ``i_interface`` is the number of the enclosed domain, starting from 0. + """ + result = np.zeros((2, 2, len(self.interfaces)), dtype=np.object) + + i_unknown = 0 + for side in self.sides: + for field_kind in self.field_kinds: + for i_interface in range(len(self.interfaces)): + + if self.is_field_present(field_kind): + dens = unknown[i_unknown] + i_unknown += 1 + else: + dens = 0 + + _, _, interface_id = self.interfaces[i_interface] + + if not with_l2_weights: + dens = sym.cse( + dens/self.get_sqrt_weight(interface_id), + "dens_{side}_{field}_{dom}".format( + side={ + self.side_out: "o", + self.side_in: "i"} + [side], + field={ + self.field_kind_e: "E", + self.field_kind_h: "H" + } + [field_kind], + dom=i_interface)) + + result[side, field_kind, i_interface] = dens + + assert i_unknown == len(unknown) + return result + + def representation(self, unknown, i_domain): + """ + :return: a symbolic expression for the representation of the PDE solution + in domain number *i_domain*. + """ + unk = self._structured_unknown(unknown, with_l2_weights=False) + + result = [] + + for field_kind in self.field_kinds: + if not self.is_field_present(field_kind): + continue + + field_result = 0 + for i_interface, (i_domain_outer, i_domain_inner, interface_id) in ( + enumerate(self.interfaces)): + if i_domain_outer == i_domain: + side = self.side_out + elif i_domain_inner == i_domain: + side = self.side_in + else: + continue + + my_unk = unk[side, field_kind, i_interface] + if my_unk: + field_result += sym.S( + self.kernel, + my_unk, + source=interface_id, + k=self.domain_K_exprs[i_domain]) + + result.append(field_result) + + from pytools.obj_array import make_obj_array + return make_obj_array(result) + + def operator(self, unknown): + density_unk = self._structured_unknown(unknown, with_l2_weights=False) + discrete_unk = self._structured_unknown(unknown, with_l2_weights=True) + + result = [] + for bc in self.bcs: + op = 0 + + for side in self.sides: + for term in bc: + unk_index = (side, term.field_kind, term.i_interface) + density = density_unk[unk_index] + discrete = discrete_unk[unk_index] + + op += self.bc_term_to_operator_contrib( + term, side, sym.S, density, discrete) + + result.append(op) + + return np.array(result, dtype=np.object) + +# }}} + + +# {{{ single + double layer representation + +class DielectricSDRep2DBoundaryOperator(Dielectric2DBoundaryOperatorBase): + pot_kind_S = 0 + pot_kind_D = 1 + pot_kinds = [pot_kind_S, pot_kind_D] + potential_ops = { + pot_kind_S: sym.S, + pot_kind_D: sym.D, + } + + def __init__(self, mode, k_vacuum, domain_k_exprs, beta, + interfaces, use_l2_weighting=None): + + super(DielectricSDRep2DBoundaryOperator, self).__init__( + mode, k_vacuum, domain_k_exprs, beta, + interfaces, use_l2_weighting=use_l2_weighting) + + side_in = self.side_in + side_out = self.side_out + + def find_normal_derivative_bc_coeff(field_kind, i_interface, side): + result = 0 + for bc in self.bcs: + for term in bc: + if term.field_kind != field_kind: + continue + if term.i_interface != i_interface: + continue + if term.direction != self.dir_normal: + continue + + if side == side_in: + result += term.coeff_inner + elif side == side_out: + result += term.coeff_outer + else: + raise ValueError("invalid side") + + return result + + self.density_coeffs = np.zeros( + (len(self.pot_kinds), len(self.field_kinds), + len(self.interfaces), len(self.sides)), + dtype=np.object) + for field_kind in self.field_kinds: + for i_interface in range(len(self.interfaces)): + self.density_coeffs[ + self.pot_kind_S, field_kind, i_interface, side_in] = 1 + self.density_coeffs[ + self.pot_kind_S, field_kind, i_interface, side_out] = 1 + + # These need to satisfy + # + # [dens_coeff_D * bc_coeff * dn D] + # = dens_coeff_D_out * bc_coeff_out * (dn D) + # + dens_coeff_D_in * bc_coeff_in * dn D + # = 0 + # + # (because dn D is hypersingular, which we'd like to cancel out) + # + # NB: bc_coeff_{in,out} already contain the signs to realize + # the subtraction for the jump. (So the "+" above is as it + # should be.) + + dens_coeff_D_in = find_normal_derivative_bc_coeff( # noqa + field_kind, i_interface, side_out) + dens_coeff_D_out = - find_normal_derivative_bc_coeff( # noqa + field_kind, i_interface, side_in) + + self.density_coeffs[ + self.pot_kind_D, field_kind, i_interface, side_in] \ + = dens_coeff_D_in + self.density_coeffs[ + self.pot_kind_D, field_kind, i_interface, side_out] \ + = dens_coeff_D_out + + def _structured_unknown(self, unknown, with_l2_weights): + """ + :arg with_l2_weights: If True, return the 'bare' unknowns + that do not have the :math:`L^2` weights divided out. + Note: Those unknowns should *not* be interpreted as + point values of a density. + :returns: an array of unknowns, with the following index axes: + ``[pot_kind, field_kind, i_interface]``, where + ``pot_kind`` is 0 for the single-layer part and 1 for the double-layer + part, + ``field_kind`` is 0 for the E-field and 1 for the H-field part, + ``i_interface`` is the number of the enclosed domain, starting from 0. + """ + result = np.zeros((2, 2, len(self.interfaces)), dtype=np.object) + + i_unknown = 0 + for pot_kind in self.pot_kinds: + for field_kind in self.field_kinds: + for i_interface in range(len(self.interfaces)): + + if self.is_field_present(field_kind): + dens = unknown[i_unknown] + i_unknown += 1 + else: + dens = 0 + + _, _, interface_id = self.interfaces[i_interface] + + if not with_l2_weights: + dens = sym.cse( + dens/self.get_sqrt_weight(interface_id), + "dens_{pot}_{field}_{intf}".format( + pot={0: "S", 1: "D"}[pot_kind], + field={ + self.field_kind_e: "E", + self.field_kind_h: "H" + } + [field_kind], + intf=i_interface)) + + result[pot_kind, field_kind, i_interface] = dens + + assert i_unknown == len(unknown) + return result + + def representation(self, unknown, i_domain): + """ + :return: a symbolic expression for the representation of the PDE solution + in domain number *i_domain*. + """ + unk = self._structured_unknown(unknown, with_l2_weights=False) + + result = [] + + for field_kind in self.field_kinds: + if not self.is_field_present(field_kind): + continue + + field_result = 0 + for pot_kind in self.pot_kinds: + for i_interface, (i_domain_outer, i_domain_inner, interface_id) in ( + enumerate(self.interfaces)): + if i_domain_outer == i_domain: + side = self.side_out + elif i_domain_inner == i_domain: + side = self.side_in + else: + continue + + my_unk = unk[pot_kind, field_kind, i_interface] + if my_unk: + field_result += ( + self.density_coeffs[ + pot_kind, field_kind, i_interface, side] + * self.potential_ops[pot_kind]( + self.kernel, + my_unk, + source=interface_id, + k=self.domain_K_exprs[i_domain] + )) + + result.append(field_result) + + from pytools.obj_array import make_obj_array + return make_obj_array(result) + + def operator(self, unknown): + density_unk = self._structured_unknown(unknown, with_l2_weights=False) + discrete_unk = self._structured_unknown(unknown, with_l2_weights=True) + + result = [] + for bc in self.bcs: + op = 0 + + for pot_kind in self.pot_kinds: + for term in bc: + + for side in self.sides: + raw_potential_op = \ + self.potential_ops[pot_kind] + + unk_index = (pot_kind, term.field_kind, term.i_interface) + density = density_unk[unk_index] + discrete = discrete_unk[unk_index] + + op += ( + self.density_coeffs[ + pot_kind, term.field_kind, term.i_interface, + side] + * self.bc_term_to_operator_contrib( + term, side, raw_potential_op, density, discrete) + ) + + result.append(op) + + return np.array(result, dtype=np.object) + +# }}} + +# }}} + diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py index ff5a22e2..25f72db9 100644 --- a/pytential/symbolic/pde/scalar.py +++ b/pytential/symbolic/pde/scalar.py @@ -289,630 +289,4 @@ class NeumannOperator(L2WeightedPDEOperator): # }}} -# {{{ dielectric - -class Dielectric2DBoundaryOperatorBase(L2WeightedPDEOperator): - r""" - Solves the following system of BVPs on :math:`\mathbb{R}^2`, in which - a disjoint family of domains :math:`\Omega_i` is embedded: - - .. math:: - - \triangle E + (k_i^2-\beta^2) E = 0\quad \text{on $\Omega_i$}\\ - \triangle H + (k_i^2-\beta^2) H = 0\quad \text{on $\Omega_i$}\\ - [H] = 0 \text{ on $\partial \Omega_i$},\quad - [E] = 0 \text{ on $\partial \Omega_i$}\\ - \left[ \frac{k_0}{k^2-\beta^2} \partial_{\hat n}H\right] = 0 - \quad\text{on $\partial \Omega_i$},\quad\\ - \left[ \frac{k_0}{k^2-\beta^2} \partial_{\hat n}E\right] = 0 - \quad\text{on $\partial \Omega_i$} - - :math:`E` and :math:`H` are assumed to be of the form - - .. math:: - - E(x,y,z,t)=E(x,y)e^{i(\beta z-\omega t) - H(x,y,z,t)=H(x,y)e^{i(\beta z-\omega t) - - where :math:`[\cdot]` denotes the jump across an interface, and :math:`k` - (without an index) denotes the value of :math:`k` on either side of the - interface, for the purpose of computing the jump. :math:`\hat n` denotes - the unit normal of the interface. - - .. automethod:: make_unknown - .. automethod:: representation_outer - .. automethod:: representation_inner - .. automethod:: operator - """ - - field_kind_e = 0 - field_kind_h = 1 - field_kinds = [field_kind_e, field_kind_h] - - side_in = 0 - side_out = 1 - sides = [side_in, side_out] - side_to_sign = { - side_in: -1, - side_out: 1, - } - - dir_none = 0 - dir_normal = 1 - dir_tangential = 2 - - BCTermDescriptor = namedtuple("BCDescriptor", - "i_interface direction field_kind coeff_inner coeff_outer".split()) - - # {{{ constructor - - def __init__(self, mode, k_vacuum, domain_k_exprs, beta, - interfaces, use_l2_weighting=None): - """ - :attr mode: one of 'te', 'tm', 'tem' - :attr k_vacuum: A symbolic expression for the wave number in vacuum. - May be a string, which will be interpreted as a variable name. - :attr interfaces: a tuple of tuples - ``(outer_domain, inner_domain, interface_id)``, - where *outer_domain* and *inner_domain* are indices into - *domain_k_names*, - and *interface_id* is a symbolic name for the discretization of the - interface. 'outer' designates the side of the interface to which - the normal points. - :attr domain_k_exprs: a tuple of variable names of the Helmholtz - parameter *k*, to be used inside each part of the source geometry. - May also be a tuple of strings, which will be transformed into - variable references of the corresponding names. - :attr beta: A symbolic expression for the wave number in the :math:`z` - direction. May be a string, which will be interpreted as a variable - name. - """ - - if use_l2_weighting is None: - use_l2_weighting = False - - super(Dielectric2DBoundaryOperatorBase, self).__init__( - use_l2_weighting=use_l2_weighting) - - if mode == "te": - self.ez_enabled = False - self.hz_enabled = True - elif mode == "tm": - self.ez_enabled = True - self.hz_enabled = False - elif mode == "tem": - self.ez_enabled = True - self.hz_enabled = True - else: - raise ValueError("invalid mode '%s'" % mode) - - self.interfaces = interfaces - - fk_e = self.field_kind_e - fk_h = self.field_kind_h - - dir_none = self.dir_none - dir_normal = self.dir_normal - dir_tangential = self.dir_tangential - - if isinstance(beta, str): - beta = sym.var(beta) - beta = sym.cse(beta, "beta") - - if isinstance(k_vacuum, str): - k_vacuum = sym.var(k_vacuum) - k_vacuum = sym.cse(k_vacuum, "k_vac") - - self.domain_k_exprs = [ - sym.var(k_expr) - if isinstance(k_expr, str) - else sym.cse(k_expr, "k%d" % idom) - for idom, k_expr in enumerate(domain_k_exprs)] - del domain_k_exprs - - # Note the case of k/K! - # "K" is the 2D Helmholtz parameter. - # "k" is the 3D Helmholtz parameter. - - self.domain_K_exprs = [ - sym.cse((k_expr**2-beta**2)**0.5, "K%d" % i) - for i, k_expr in enumerate(self.domain_k_exprs)] - - from sumpy.kernel import HelmholtzKernel - self.kernel = HelmholtzKernel(2, allow_evanescent=True) - - # {{{ build bc list - - # list of tuples, where each tuple consists of BCTermDescriptor instances - - all_bcs = [] - for i_interface, (outer_domain, inner_domain, _) in ( - enumerate(self.interfaces)): - k_outer = self.domain_k_exprs[outer_domain] - k_inner = self.domain_k_exprs[inner_domain] - - all_bcs += [ - ( # [E] = 0 - self.BCTermDescriptor( - i_interface=i_interface, - direction=dir_none, - field_kind=fk_e, - coeff_outer=1, - coeff_inner=-1), - ), - ( # [H] = 0 - self.BCTermDescriptor( - i_interface=i_interface, - direction=dir_none, - field_kind=fk_h, - coeff_outer=1, - coeff_inner=-1), - ), - ( - self.BCTermDescriptor( - i_interface=i_interface, - direction=dir_tangential, - field_kind=fk_e, - coeff_outer=beta/(k_outer**2-beta**2), - coeff_inner=-beta/(k_inner**2-beta**2)), - self.BCTermDescriptor( - i_interface=i_interface, - direction=dir_normal, - field_kind=fk_h, - coeff_outer=sym.cse(-k_vacuum/(k_outer**2-beta**2)), - coeff_inner=sym.cse(k_vacuum/(k_inner**2-beta**2))), - ), - ( - self.BCTermDescriptor( - i_interface=i_interface, - direction=dir_tangential, - field_kind=fk_h, - coeff_outer=beta/(k_outer**2-beta**2), - coeff_inner=-beta/(k_inner**2-beta**2)), - self.BCTermDescriptor( - i_interface=i_interface, - direction=dir_normal, - field_kind=fk_e, - coeff_outer=sym.cse( - (k_outer**2/k_vacuum)/(k_outer**2-beta**2)), - coeff_inner=sym.cse( - -(k_inner**2/k_vacuum) - / (k_inner**2-beta**2))) - ), - ] - - del k_outer - del k_inner - - self.bcs = [] - for bc in all_bcs: - any_significant_e = any( - term.field_kind == fk_e - and term.direction in [dir_normal, dir_none] - for term in bc) - any_significant_h = any( - term.field_kind == fk_h - and term.direction in [dir_normal, dir_none] - for term in bc) - is_necessary = ( - (self.ez_enabled and any_significant_e) - or - (self.hz_enabled and any_significant_h)) - - # Only keep tangential modes for TEM. Otherwise, - # no jump in H already implies jump condition on - # tangential derivative. - is_tem = self.ez_enabled and self.hz_enabled - terms = tuple( - term - for term in bc - if term.direction != dir_tangential - or is_tem) - - if is_necessary: - self.bcs.append(terms) - - assert (len(all_bcs) - * (int(self.ez_enabled) + int(self.hz_enabled)) // 2 - == len(self.bcs)) - - # }}} - - # }}} - - def is_field_present(self, field_kind): - return ( - (field_kind == self.field_kind_e and self.ez_enabled) - or - (field_kind == self.field_kind_h and self.hz_enabled)) - - def make_unknown(self, name): - num_densities = ( - 2 - * (int(self.ez_enabled) + int(self.hz_enabled)) - * len(self.interfaces)) - - assert num_densities == len(self.bcs) - - return sym.make_sym_vector(name, num_densities) - - def bc_term_to_operator_contrib(self, term, side, raw_potential_op, - density, discrete): - potential_op = raw_potential_op - - side_sign = self.side_to_sign[side] - - domain_outer, domain_inner, interface_id = \ - self.interfaces[term.i_interface] - if side == self.side_in: - K_expr = self.domain_K_exprs[domain_inner] # noqa - bc_coeff = term.coeff_inner - elif side == self.side_out: - K_expr = self.domain_K_exprs[domain_outer] # noqa - bc_coeff = term.coeff_outer - else: - raise ValueError("invalid value of 'side'") - - potential_op = potential_op( - self.kernel, density, source=interface_id, - k=K_expr) - - if term.direction == self.dir_none: - if raw_potential_op is sym.S: - jump_term = 0 - elif raw_potential_op is sym.D: - jump_term = (side_sign*0.5) * discrete - else: - assert False, raw_potential_op - elif term.direction == self.dir_normal: - potential_op = sym.normal_derivative( - potential_op, interface_id) - - if raw_potential_op is sym.S: - # S' - jump_term = (-side_sign*0.5) * discrete - elif raw_potential_op is sym.D: - jump_term = 0 - else: - assert False, raw_potential_op - - elif term.direction == self.dir_tangential: - potential_op = sym.tangential_derivative( - raw_potential_op( - self.kernel, density, source=interface_id, - k=K_expr, qbx_forced_limit=side_sign), - interface_id).a.as_scalar() - - # Some of these may have jumps, but QBX does the dirty - # work here by directly computing the limit. - jump_term = 0 - - else: - raise ValueError("invalid direction") - - potential_op = ( - jump_term - + self.get_sqrt_weight(interface_id)*potential_op) - - del jump_term - - contrib = bc_coeff * potential_op - - if (raw_potential_op is sym.D - and term.direction == self.dir_normal): - # FIXME The hypersingular part should perhaps be - # treated specially to avoid cancellation. - pass - - return contrib - - -# {{{ single-layer representation - -class DielectricSRep2DBoundaryOperator(Dielectric2DBoundaryOperatorBase): - def _structured_unknown(self, unknown, with_l2_weights): - """ - :arg with_l2_weights: If True, return the 'bare' unknowns - that do not have the :math:`L^2` weights divided out. - Note: Those unknowns should *not* be interpreted as - point values of a density. - :returns: an array of unknowns, with the following index axes: - ``[side, field_kind, i_interface]``, where - ``side`` is o for the outside part and i for the interior part, - ``field_kind`` is 0 for the E-field and 1 for the H-field part, - ``i_interface`` is the number of the enclosed domain, starting from 0. - """ - result = np.zeros((2, 2, len(self.interfaces)), dtype=np.object) - - i_unknown = 0 - for side in self.sides: - for field_kind in self.field_kinds: - for i_interface in range(len(self.interfaces)): - - if self.is_field_present(field_kind): - dens = unknown[i_unknown] - i_unknown += 1 - else: - dens = 0 - - _, _, interface_id = self.interfaces[i_interface] - - if not with_l2_weights: - dens = sym.cse( - dens/self.get_sqrt_weight(interface_id), - "dens_{side}_{field}_{dom}".format( - side={ - self.side_out: "o", - self.side_in: "i"} - [side], - field={ - self.field_kind_e: "E", - self.field_kind_h: "H" - } - [field_kind], - dom=i_interface)) - - result[side, field_kind, i_interface] = dens - - assert i_unknown == len(unknown) - return result - - def representation(self, unknown, i_domain): - """ - :return: a symbolic expression for the representation of the PDE solution - in domain number *i_domain*. - """ - unk = self._structured_unknown(unknown, with_l2_weights=False) - - result = [] - - for field_kind in self.field_kinds: - if not self.is_field_present(field_kind): - continue - - field_result = 0 - for i_interface, (i_domain_outer, i_domain_inner, interface_id) in ( - enumerate(self.interfaces)): - if i_domain_outer == i_domain: - side = self.side_out - elif i_domain_inner == i_domain: - side = self.side_in - else: - continue - - my_unk = unk[side, field_kind, i_interface] - if my_unk: - field_result += sym.S( - self.kernel, - my_unk, - source=interface_id, - k=self.domain_K_exprs[i_domain]) - - result.append(field_result) - - from pytools.obj_array import make_obj_array - return make_obj_array(result) - - def operator(self, unknown): - density_unk = self._structured_unknown(unknown, with_l2_weights=False) - discrete_unk = self._structured_unknown(unknown, with_l2_weights=True) - - result = [] - for bc in self.bcs: - op = 0 - - for side in self.sides: - for term in bc: - unk_index = (side, term.field_kind, term.i_interface) - density = density_unk[unk_index] - discrete = discrete_unk[unk_index] - - op += self.bc_term_to_operator_contrib( - term, side, sym.S, density, discrete) - - result.append(op) - - return np.array(result, dtype=np.object) - -# }}} - - -# {{{ single + double layer representation - -class DielectricSDRep2DBoundaryOperator(Dielectric2DBoundaryOperatorBase): - pot_kind_S = 0 - pot_kind_D = 1 - pot_kinds = [pot_kind_S, pot_kind_D] - potential_ops = { - pot_kind_S: sym.S, - pot_kind_D: sym.D, - } - - def __init__(self, mode, k_vacuum, domain_k_exprs, beta, - interfaces, use_l2_weighting=None): - - super(DielectricSDRep2DBoundaryOperator, self).__init__( - mode, k_vacuum, domain_k_exprs, beta, - interfaces, use_l2_weighting=use_l2_weighting) - - side_in = self.side_in - side_out = self.side_out - - def find_normal_derivative_bc_coeff(field_kind, i_interface, side): - result = 0 - for bc in self.bcs: - for term in bc: - if term.field_kind != field_kind: - continue - if term.i_interface != i_interface: - continue - if term.direction != self.dir_normal: - continue - - if side == side_in: - result += term.coeff_inner - elif side == side_out: - result += term.coeff_outer - else: - raise ValueError("invalid side") - - return result - - self.density_coeffs = np.zeros( - (len(self.pot_kinds), len(self.field_kinds), - len(self.interfaces), len(self.sides)), - dtype=np.object) - for field_kind in self.field_kinds: - for i_interface in range(len(self.interfaces)): - self.density_coeffs[ - self.pot_kind_S, field_kind, i_interface, side_in] = 1 - self.density_coeffs[ - self.pot_kind_S, field_kind, i_interface, side_out] = 1 - - # These need to satisfy - # - # [dens_coeff_D * bc_coeff * dn D] - # = dens_coeff_D_out * bc_coeff_out * (dn D) - # + dens_coeff_D_in * bc_coeff_in * dn D - # = 0 - # - # (because dn D is hypersingular, which we'd like to cancel out) - # - # NB: bc_coeff_{in,out} already contain the signs to realize - # the subtraction for the jump. (So the "+" above is as it - # should be.) - - dens_coeff_D_in = find_normal_derivative_bc_coeff( # noqa - field_kind, i_interface, side_out) - dens_coeff_D_out = - find_normal_derivative_bc_coeff( # noqa - field_kind, i_interface, side_in) - - self.density_coeffs[ - self.pot_kind_D, field_kind, i_interface, side_in] \ - = dens_coeff_D_in - self.density_coeffs[ - self.pot_kind_D, field_kind, i_interface, side_out] \ - = dens_coeff_D_out - - def _structured_unknown(self, unknown, with_l2_weights): - """ - :arg with_l2_weights: If True, return the 'bare' unknowns - that do not have the :math:`L^2` weights divided out. - Note: Those unknowns should *not* be interpreted as - point values of a density. - :returns: an array of unknowns, with the following index axes: - ``[pot_kind, field_kind, i_interface]``, where - ``pot_kind`` is 0 for the single-layer part and 1 for the double-layer - part, - ``field_kind`` is 0 for the E-field and 1 for the H-field part, - ``i_interface`` is the number of the enclosed domain, starting from 0. - """ - result = np.zeros((2, 2, len(self.interfaces)), dtype=np.object) - - i_unknown = 0 - for pot_kind in self.pot_kinds: - for field_kind in self.field_kinds: - for i_interface in range(len(self.interfaces)): - - if self.is_field_present(field_kind): - dens = unknown[i_unknown] - i_unknown += 1 - else: - dens = 0 - - _, _, interface_id = self.interfaces[i_interface] - - if not with_l2_weights: - dens = sym.cse( - dens/self.get_sqrt_weight(interface_id), - "dens_{pot}_{field}_{intf}".format( - pot={0: "S", 1: "D"}[pot_kind], - field={ - self.field_kind_e: "E", - self.field_kind_h: "H" - } - [field_kind], - intf=i_interface)) - - result[pot_kind, field_kind, i_interface] = dens - - assert i_unknown == len(unknown) - return result - - def representation(self, unknown, i_domain): - """ - :return: a symbolic expression for the representation of the PDE solution - in domain number *i_domain*. - """ - unk = self._structured_unknown(unknown, with_l2_weights=False) - - result = [] - - for field_kind in self.field_kinds: - if not self.is_field_present(field_kind): - continue - - field_result = 0 - for pot_kind in self.pot_kinds: - for i_interface, (i_domain_outer, i_domain_inner, interface_id) in ( - enumerate(self.interfaces)): - if i_domain_outer == i_domain: - side = self.side_out - elif i_domain_inner == i_domain: - side = self.side_in - else: - continue - - my_unk = unk[pot_kind, field_kind, i_interface] - if my_unk: - field_result += ( - self.density_coeffs[ - pot_kind, field_kind, i_interface, side] - * self.potential_ops[pot_kind]( - self.kernel, - my_unk, - source=interface_id, - k=self.domain_K_exprs[i_domain] - )) - - result.append(field_result) - - from pytools.obj_array import make_obj_array - return make_obj_array(result) - - def operator(self, unknown): - density_unk = self._structured_unknown(unknown, with_l2_weights=False) - discrete_unk = self._structured_unknown(unknown, with_l2_weights=True) - - result = [] - for bc in self.bcs: - op = 0 - - for pot_kind in self.pot_kinds: - for term in bc: - - for side in self.sides: - raw_potential_op = \ - self.potential_ops[pot_kind] - - unk_index = (pot_kind, term.field_kind, term.i_interface) - density = density_unk[unk_index] - discrete = discrete_unk[unk_index] - - op += ( - self.density_coeffs[ - pot_kind, term.field_kind, term.i_interface, - side] - * self.bc_term_to_operator_contrib( - term, side, raw_potential_op, density, discrete) - ) - - result.append(op) - - return np.array(result, dtype=np.object) - -# }}} - -# }}} - # vim: foldmethod=marker -- GitLab From e175468e5d3c0818a5a13fad07dbead195384630 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 13 Dec 2016 17:53:11 -0600 Subject: [PATCH 2/8] Source derivative direction vector finding: support quotients --- pytential/symbolic/primitives.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 5e8ff991..2e56c83d 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -585,6 +585,9 @@ def _get_dir_vec(dsource, ambient_dim): def map_common_subexpression(self, expr): return {1: expr} + def map_quotient(self, expr): + return {1: expr} + coeffs = _DSourceCoefficientFinder()(dsource) dir_vec = np.zeros(ambient_dim, np.object) -- GitLab From 09af905acba932243ca985e48a43b295654a2574 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 13 Dec 2016 17:54:14 -0600 Subject: [PATCH 3/8] Solve: clean up after python-modernize --- pytential/solve.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pytential/solve.py b/pytential/solve.py index fa993642..34081d27 100644 --- a/pytential/solve.py +++ b/pytential/solve.py @@ -1,7 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -from __future__ import print_function -from six.moves import range +from __future__ import division, absolute_import, print_function __copyright__ = "Copyright (C) 2012-2013 Andreas Kloeckner" @@ -26,6 +23,8 @@ THE SOFTWARE. """ +from six.moves import range + __doc__ = """ .. autofunction:: gmres -- GitLab From 7e3b334ace0ba5a092408c46d82b2626856755f6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 13 Dec 2016 17:54:43 -0600 Subject: [PATCH 4/8] Mode finding SKIE: Operator coded up [ci skip] --- examples/find-photonic-mode-sk.py | 55 +++++++++++------ pytential/symbolic/pde/maxwell/fiber.py | 81 ++++++++++++------------- 2 files changed, 75 insertions(+), 61 deletions(-) diff --git a/examples/find-photonic-mode-sk.py b/examples/find-photonic-mode-sk.py index c44a6ce8..9703b868 100644 --- a/examples/find-photonic-mode-sk.py +++ b/examples/find-photonic-mode-sk.py @@ -25,17 +25,20 @@ THE SOFTWARE. """ import numpy as np -import numpy.linalg as la import pyopencl as cl -from pytential import sym +from pytential import sym, bind from functools import partial from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) +import logging + def find_mode(): + logging.basicConfig(level=logging.INFO) + import warnings warnings.simplefilter("error", np.ComplexWarning) @@ -44,21 +47,25 @@ def find_mode(): k0 = 1.4447 k1 = k0*1.02 - beta_sym = sym.var("beta") from pytential.symbolic.pde.maxwell.fiber import \ SecondKindInfZMuellerOperator pde_op = SecondKindInfZMuellerOperator( - k_vacuum=1, + k_vacuum="k_v", interfaces=((0, 1, sym.DEFAULT_SOURCE),), - domain_k_exprs=(k0, k1), - beta=beta_sym, + domain_k_exprs=("k0", "k1"), + beta="beta", use_l2_weighting=True) - jm_sym = pde_op.make_unknown("") - op = pde_op.operator(jm_sym) - 1/0 + base_context = { + "k0": k0, + "k1": k1, + "k_v": 1, + } + + u_sym = pde_op.make_unknown("u") + op = pde_op.operator(u_sym) # {{{ discretization setup @@ -66,7 +73,7 @@ def find_mode(): curve_f = partial(ellipse, 1) target_order = 7 - qbx_order = 4 + qbx_order = 3 nelements = 30 from meshmode.mesh.processing import affine_map @@ -81,27 +88,37 @@ def find_mode(): from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory from pytential.qbx import QBXLayerPotentialSource - density_discr = Discretization( + pre_density_discr = Discretization( cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) - qbx = QBXLayerPotentialSource(density_discr, 4*target_order, + qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, qbx_order, # Don't use FMM for now - fmm_order=False) + fmm_order=qbx_order+5).with_refinement() + density_discr = qbx.density_discr # }}} x_vec = np.random.randn(len(u_sym)*density_discr.nnodes) y_vec = np.random.randn(len(u_sym)*density_discr.nnodes) - def muller_solve_func(beta): - from pytential.symbolic.execution import build_matrix - mat = build_matrix( - queue, qbx, op, u_sym, - context={"beta": beta}).get() + bound_op = bind(qbx, op) - return 1/x_vec.dot(la.solve(mat, y_vec)) + def muller_solve_func(beta): + from pytential.solve import gmres + gmres_result = gmres( + bound_op.scipy_op(queue, "u", + np.complex128, beta=beta, **base_context), + y_vec, tol=1e-12, progress=True, + stall_iterations=0, + hard_failure=False) + minv_y = gmres_result.solution + print("gmres state:", gmres_result.state) + + z = 1/x_vec.dot(minv_y) + print("muller func value:", repr(z)) + return z starting_guesses = (1+0j)*( k0 diff --git a/pytential/symbolic/pde/maxwell/fiber.py b/pytential/symbolic/pde/maxwell/fiber.py index cd11ca63..6b746bdf 100644 --- a/pytential/symbolic/pde/maxwell/fiber.py +++ b/pytential/symbolic/pde/maxwell/fiber.py @@ -41,9 +41,6 @@ from collections import namedtuple from six.moves import range from pytential import sym -from pytential.symbolic.primitives import ( - cse, - sqrt_jac_q_weight, QWeight, area_element) from pytools import memoize_method from pytential.symbolic.pde.scalar import L2WeightedPDEOperator @@ -79,18 +76,14 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): self.interfaces = interfaces - if isinstance(beta, str): - beta = sym.var(beta) + beta = sym.var(beta) self.beta = sym.cse(beta, "beta") - if isinstance(k_vacuum, str): - k_vacuum = sym.var(k_vacuum) + k_vacuum = sym.var(k_vacuum) self.k_vacuum = sym.cse(k_vacuum, "k_vac") self.domain_k_exprs = [ sym.var(k_expr) - if isinstance(k_expr, str) - else sym.cse(k_expr, "k%d" % idom) for idom, k_expr in enumerate(domain_k_exprs)] del domain_k_exprs @@ -101,17 +94,8 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): from sumpy.kernel import HelmholtzKernel self.kernel = HelmholtzKernel(2, allow_evanescent=True) - def make_unknown(self, prefix): - from pytools import reverse_dictionary - index_to_name_and_intf = reverse_dictionary(self.unknown_index) - - from pytools.obj_array import make_obj_array - return make_obj_array([ - sym.var(prefix - + index_to_name_and_intf[i][0] - + str(index_to_name_and_intf[i][1])) - for i in range(len(index_to_name_and_intf)) - ]) + def make_unknown(self, name): + return sym.make_sym_vector(name, len(self.unknown_index)) @property @memoize_method @@ -128,7 +112,7 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): raise NotImplementedError() def operator(self, unknown): - result = np.zeros(4*len(self.interfaces)) + result = np.zeros(4*len(self.interfaces), dtype=object) unk_idx = self.unknown_index @@ -138,47 +122,56 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): idx_mt = unk_idx["mt", i] idx_mz = unk_idx["mz", i] - phi1 = jt = unknown[idx_jt] - phi2 = jz = unknown[idx_jz] - phi3 = mt = unknown[idx_mt] - phi4 = mz = unknown[idx_mz] + phi1 = unknown[idx_jt] + phi2 = unknown[idx_jz] + phi3 = unknown[idx_mt] + phi4 = unknown[idx_mz] ne = self.beta/self.k_vacuum dom0_idx, dom1_idx, where = self.interfaces[i] dom_indices = [dom0_idx, dom1_idx] - def S(dom, density): # noqa + tangent = (sym.pseudoscalar(2, 1, where) + / sym.area_element(2, 1, where)) + normal = sym.normal(2, 1, where) + + def S(dom, density, qbx_forced_limit=+1): # noqa return sym.S( self.kernel, density, - k=self.domain_K_exprs[dom_indices[dom]]) + k=self.domain_K_exprs[dom_indices[dom]], + qbx_forced_limit=qbx_forced_limit) - def D(dom, density): # noqa + def D(dom, density, qbx_forced_limit="avg"): # noqa return sym.D( self.kernel, density, - k=self.domain_K_exprs[dom_indices[dom]]) + k=self.domain_K_exprs[dom_indices[dom]], + qbx_forced_limit=qbx_forced_limit) def T(dom, density): # noqa - 1/0 + return sym.int_g_dsource( + 2, + tangent, + self.kernel, + density, + k=self.domain_K_exprs[dom_indices[dom]], + # ???? + qbx_forced_limit="avg").xproject(0) def Tt(dom, density): # noqa - 1/0 + return sym.tangential_derivative(2, T(dom, density)).xproject(0) def Sn(dom, density): # noqa - 1/0 + return sym.normal_derivative( + 2, + S(dom, density, + qbx_forced_limit="avg")) def St(dom, density): # noqa - 1/0 + return sym.tangential_derivative(2, S(dom, density)).xproject(0) - tangent = (sym.pseudoscalar(2, 1, where) - / sym.area_element(2, 1, where)) - normal = sym.normal(2, 1, where) - print(sym.pretty( - )) - 1/0 - - n0 = 1 # FIXME - n1 = 2 # FIXME + n0 = self.domain_k_exprs[dom0_idx] # FIXME + n1 = self.domain_k_exprs[dom1_idx] # FIXME a11 = n0**2 * D(0, phi1) - n1**2 * D(1, phi1) a22 = -n0**2 * Sn(0, phi2) + n1**2 * Sn(1, phi2) @@ -225,6 +218,9 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): result[idx_mz] += d4 + a41 + a42 + a43 + a44 # TODO: L2 weighting + # TODO: Add representation contributions to other boundaries + # abutting the domain + return result # }}} @@ -855,3 +851,4 @@ class DielectricSDRep2DBoundaryOperator(Dielectric2DBoundaryOperatorBase): # }}} +# vim: foldmethod=marker -- GitLab From aae65cc9a686335b44a79bd85d298acd81d60d7f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 14 Dec 2016 01:03:36 -0600 Subject: [PATCH 5/8] Waveguide SKIE: Minor fixes --- pytential/symbolic/pde/maxwell/fiber.py | 40 ++++++++++++++----------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/fiber.py b/pytential/symbolic/pde/maxwell/fiber.py index 6b746bdf..0f42ed8a 100644 --- a/pytential/symbolic/pde/maxwell/fiber.py +++ b/pytential/symbolic/pde/maxwell/fiber.py @@ -117,24 +117,28 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): unk_idx = self.unknown_index for i in range(len(self.interfaces)): - idx_jt = unk_idx["jt", i] - idx_jz = unk_idx["jz", i] - idx_mt = unk_idx["mt", i] - idx_mz = unk_idx["mz", i] + idx_jt = unk_idx["jz", i] + idx_jz = unk_idx["jt", i] + idx_mt = unk_idx["mz", i] + idx_mz = unk_idx["mt", i] phi1 = unknown[idx_jt] phi2 = unknown[idx_jz] phi3 = unknown[idx_mt] phi4 = unknown[idx_mz] - ne = self.beta/self.k_vacuum + ne = sym.cse(self.beta/self.k_vacuum, "ne") dom0_idx, dom1_idx, where = self.interfaces[i] dom_indices = [dom0_idx, dom1_idx] - tangent = (sym.pseudoscalar(2, 1, where) - / sym.area_element(2, 1, where)) - normal = sym.normal(2, 1, where) + tangent = sym.cse( + (sym.pseudoscalar(2, 1, where) + / sym.area_element(2, 1, where)), + "tangent") + normal = sym.cse( + sym.normal(2, 1, where), + "normal") def S(dom, density, qbx_forced_limit=+1): # noqa return sym.S( @@ -156,7 +160,7 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): density, k=self.domain_K_exprs[dom_indices[dom]], # ???? - qbx_forced_limit="avg").xproject(0) + qbx_forced_limit=1).xproject(0) def Tt(dom, density): # noqa return sym.tangential_derivative(2, T(dom, density)).xproject(0) @@ -170,13 +174,15 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): def St(dom, density): # noqa return sym.tangential_derivative(2, S(dom, density)).xproject(0) - n0 = self.domain_k_exprs[dom0_idx] # FIXME - n1 = self.domain_k_exprs[dom1_idx] # FIXME + n0 = sym.cse(self.domain_k_exprs[dom0_idx]/self.k_vacuum, + "n%d" % dom0_idx) + n1 = sym.cse(self.domain_k_exprs[dom1_idx]/self.k_vacuum, + "n%d" % dom1_idx) - a11 = n0**2 * D(0, phi1) - n1**2 * D(1, phi1) - a22 = -n0**2 * Sn(0, phi2) + n1**2 * Sn(1, phi2) - a33 = D(0, phi3)-D(1, phi3) - a44 = -Sn(0, phi4) + Sn(1, phi4) + a11 = sym.cse(n0**2 * D(0, phi1) - n1**2 * D(1, phi1), "a11") + a22 = sym.cse(-n0**2 * Sn(0, phi2) + n1**2 * Sn(1, phi2), "a22") + a33 = sym.cse(D(0, phi3)-D(1, phi3), "a33") + a44 = sym.cse(-Sn(0, phi4) + Sn(1, phi4), "a44") a21 = -1j * ne * ( n0**2 * tangent.scalar_product( @@ -207,8 +213,8 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): a41 = -a23 a42 = -a24 - d1 = (n0**2 + n1**2)*phi1 - d2 = (n0**2 + n1**2)*phi2 + d1 = (n0**2 + n1**2)/2 * phi1 + d2 = (n0**2 + n1**2)/2 * phi2 d3 = phi3 d4 = phi4 -- GitLab From 8e30ce732c3dde00b8e67e1bf7b218b1d3f4301a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 14 Dec 2016 09:23:35 -0600 Subject: [PATCH 6/8] Waveguide tweaks --- examples/find-photonic-mode-sk.py | 11 +++++--- .../pde/maxwell/{fiber.py => waveguide.py} | 27 ++++++++++--------- 2 files changed, 21 insertions(+), 17 deletions(-) rename pytential/symbolic/pde/maxwell/{fiber.py => waveguide.py} (98%) diff --git a/examples/find-photonic-mode-sk.py b/examples/find-photonic-mode-sk.py index 9703b868..84eda74d 100644 --- a/examples/find-photonic-mode-sk.py +++ b/examples/find-photonic-mode-sk.py @@ -48,7 +48,7 @@ def find_mode(): k0 = 1.4447 k1 = k0*1.02 - from pytential.symbolic.pde.maxwell.fiber import \ + from pytential.symbolic.pde.maxwell.waveguide import \ SecondKindInfZMuellerOperator pde_op = SecondKindInfZMuellerOperator( @@ -92,10 +92,13 @@ def find_mode(): cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) - qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, + qbx, _ = QBXLayerPotentialSource( + pre_density_discr, + 4*target_order, qbx_order, - # Don't use FMM for now - fmm_order=qbx_order+5).with_refinement() + #fmm_order=qbx_order+5 + fmm_order=False + ).with_refinement() density_discr = qbx.density_discr # }}} diff --git a/pytential/symbolic/pde/maxwell/fiber.py b/pytential/symbolic/pde/maxwell/waveguide.py similarity index 98% rename from pytential/symbolic/pde/maxwell/fiber.py rename to pytential/symbolic/pde/maxwell/waveguide.py index 0f42ed8a..eed54fbe 100644 --- a/pytential/symbolic/pde/maxwell/fiber.py +++ b/pytential/symbolic/pde/maxwell/waveguide.py @@ -24,8 +24,8 @@ THE SOFTWARE. __doc__ = """ -Second-Kind Fiberoptic -^^^^^^^^^^^^^^^^^^^^^^ +Second-Kind Waveguide +^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: SecondKindInfZMuellerOperator @@ -184,30 +184,31 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): a33 = sym.cse(D(0, phi3)-D(1, phi3), "a33") a44 = sym.cse(-Sn(0, phi4) + Sn(1, phi4), "a44") - a21 = -1j * ne * ( + a21 = sym.cse(-1j * ne * ( n0**2 * tangent.scalar_product( S(0, normal * phi1)) - n1**2 * tangent.scalar_product( - S(1, normal * phi1))) - a43 = -1j * ne * ( + S(1, normal * phi1))), "a21") + + a43 = sym.cse(-1j * ne * ( tangent.scalar_product( S(0, normal * phi3)) - tangent.scalar_product( - S(1, normal * phi3))) - a13 = ne*(T(0, phi3) - T(1, phi3)) - a24 = ne*(St(0, phi4) - St(1, phi4)) - a14 = 1j*( + S(1, normal * phi3))), "a43") + a13 = sym.cse(ne*(T(0, phi3) - T(1, phi3)), "a13") + a24 = sym.cse(ne*(St(0, phi4) - St(1, phi4)), "a24") + a14 = sym.cse(1j*( (n0**2 - ne**2) * S(0, phi4) - (n1**2 - ne**2) * S(1, phi4) - ) + ), "a14") - a23 = ( + a23 = sym.cse(( 1j * (Tt(0, phi3) - Tt(1, phi3)) - 1j * ne * ( n0**2 * tangent.scalar_product( S(0, normal * phi3)) - n1**2 * tangent.scalar_product( - S(1, normal * phi3)))) + S(1, normal * phi3)))), "a23") a31 = -a13 a32 = -a14 a41 = -a23 @@ -231,7 +232,7 @@ class SecondKindInfZMuellerOperator(L2WeightedPDEOperator): # }}} -# {{{ old-style fiber +# {{{ old-style waveguide class Dielectric2DBoundaryOperatorBase(L2WeightedPDEOperator): r""" -- GitLab From 6a460df770384d9309f9bd6ccc3a4bf3d5af05f0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 5 Jan 2017 15:16:27 +0800 Subject: [PATCH 7/8] Track maxwell rename in flake8 ignore --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 377c7dc2..3662361d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -2,6 +2,6 @@ ignore = E126,E127,E128,E123,E226,E241,E242,E265,E402,W503 max-line-length=85 exclude= - pytential/symbolic/pde/maxwell.py, + pytential/symbolic/pde/maxwell/__init__.py, pytential/symbolic/old_diffop_primitives.py, -- GitLab From 6898b61f1b034b417b74273820ea3e814332dfe6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 5 Jan 2017 15:25:33 +0800 Subject: [PATCH 8/8] Flake8 fixes --- pytential/symbolic/pde/scalar.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py index 25f72db9..15f542b7 100644 --- a/pytential/symbolic/pde/scalar.py +++ b/pytential/symbolic/pde/scalar.py @@ -40,9 +40,7 @@ from pytential import sym from pytential.symbolic.primitives import ( cse, sqrt_jac_q_weight, QWeight, area_element) -import numpy as np -from collections import namedtuple -from six.moves import range +import numpy as np # noqa # {{{ L^2 weighting -- GitLab