From 6a46fcd19b920bf93d98c2c60b7fd383300aeb9a Mon Sep 17 00:00:00 2001 From: Manas Rachh Date: Tue, 6 Jun 2017 16:18:10 -0400 Subject: [PATCH 01/11] include aug muller rep --- pytential/symbolic/pde/maxwell/__init__.py | 77 +++++++++++++++++++++- 1 file changed, 76 insertions(+), 1 deletion(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 9096998a..e9301835 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -28,7 +28,7 @@ from pymbolic.primitives import Variable from pytential import sym -# {{{ MFIE +# {{{ Charge-Current MFIE class PECAugmentedMFIEOperator: """Magnetic Field Integral Equation operator, @@ -87,6 +87,81 @@ class PECAugmentedMFIEOperator: # }}} +# {{{ Charge-Current Mueller MFIE + +class MuellerAugmentedMFIEOperator: + """Magnetic Field Integral Equation operator, + under the assumption of no surface charges. + """ + + def __init__(self,omg,mu0,eps0,mu1,eps1): + from sumpy.kernel import HelmholtzKernel + self.kernel = HelmholtzKernel(3) + self.omg = omg + self.mu0 = mu0 + self.eps0 = eps0 + self.mu1=mu1 + self.eps1=eps1 + self.k0 = omg*sqrt(eps0*mu0) + self.k1 = omg*sqrt(eps1*mu1) + + def make_unknown(self, name): + return sym.make_sym_vector(name, 6) + + unk_structure = namedtuple(["jt", "rho_e", "mt", "rho_m"] + def split_unknown(self, unk): + return self.unk_structure( + jt=unk[:2], + rho_e=unk[2], + mt=unk[3:5], + rho_m=unk[5]) + + def augmul_operator(self, unk): + u = self.split_unknown(unk) + + Jxyz = cse(tangential_to_xyz(u.jt), "Jxyz") + Mxyz = cse(tangential_to_xyz(u.mt), "Mxyz") + E0 = cse(1j*self.omg*self.mu0*self.eps0*S(self.k0,0,Jxyz) + + self.mu0*gradScross(self.k0,Mxyz)-grad_S(self.k0,u.rho_e,3),"E0") + H0 = cse(-1j*self.omg*self.mu0*self.eps0*S(self.k0,0,Mxyz) + + self.eps0*gradScross(self.k0,Jxyz)+grad_S(self.k0,u.rho_m,3),"H0") + E1 = cse(1j*self.omg*self.mu1*self.eps1*S(self.k1,0,Jxyz) + + self.mu1*gradScross(self.k1,Mxyz)-grad_S(self.k1,u.rho_e,3),"E1") + H1 = cse(-1j*self.omg*self.mu1*self.eps1*S(self.k1,0,Mxyz) + + self.eps1*gradScross(self.k1,Jxyz)+grad_S(self.k1,u.rho_m,3),"H1") + F1 = cse(xyz_to_tangential(n_cross(H1-H0) + + 0.5*(self.eps0+self.eps1)*Jxyz),"F1") + F2 = cse(n_dot(self.eps1*E1-self.eps0*E0) + + 0.5*(self.eps1+self.eps0)*u.rho_e,"F2") + F3 = cse(xyz_to_tangential(n_cross(E1-E0) + + 0.5*(self.mu0+self.mu1)*Mxyz),"F3") + # sign flip included + F4 = cse(-n_dot(self.mu1*H1-self.mu0*H0) + + 0.5*(self.mu1+self.mu0)*u.rho_m,"F4") + + return sym.join_fields([F1,F2,F3,F4]) + + def augmul_rhs(self, Einc_xyz,Hinc_xyz): + return sym.join_fields([xyz_to_tangential(n_cross(Hinc_xyz)), + n_dot(self.eps1*Einc_xyz),xyz_to_tangential(n_cross(Einc_xyz)), + n_dot(-self.mu1*Hinc_xyz)] + + def scattered_volume_field(self, sol): + u = self.split_unknown(sol) + Jxyz = sym.cse(sym.tangential_to_xyz(u.jt), "Jxyz") + Mxyz = sym.cse(sym.tangential_to_xyz(u.mt), "Mxyz") + + E0 = cse(1j*self.omg*self.mu0*self.eps0*S(self.k0,0,Jxyz) + + self.mu0*gradScross(self.k0,Mxyz)-grad_S(self.k0,u.rho_e,3),"E0") + H0 = cse(-1j*self.omg*self.mu0*self.eps0*S(self.k0,0,Mxyz) + + self.eps0*gradScross(self.k0,Jxyz)+grad_S(self.k0,u.rho_m,3),"H0") + + from pytools.obj_array import join_fields + return join_fields(E0, H0) + +# }}} + + # {{{ generalized Debye def _debye_S0_surf_div(nxnxE): -- GitLab From d96fef1183e94feb5eb29e4f845f48a74df1ed08 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 19:49:26 -0400 Subject: [PATCH 02/11] Add back some more symbolic vector calculus notation --- pytential/symbolic/primitives.py | 51 +++++++++++++++++++++++++++++--- 1 file changed, 47 insertions(+), 4 deletions(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index c4c63b99..79a48bdb 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -332,7 +332,7 @@ def normal(ambient_dim, dim=None, where=None): """Exterior unit normals.""" # Don't be tempted to add a sign here. As it is, it produces - # exterior normals for positively oriented curves. + # exterior normals for positively oriented curves and surfaces. pder = ( pseudoscalar(ambient_dim, dim, where) @@ -911,8 +911,9 @@ def tangential_onb(ambient_dim, dim=None, where=None): q = avec for j in range(k): q = q - np.dot(avec, orth_pd_mat[:, j])*orth_pd_mat[:, j] + q = cse(q, "q%d" % k) - orth_pd_mat[:, k] = cse(q/sqrt(np.sum(q**2)), "orth_pd_vec") + orth_pd_mat[:, k] = cse(q/sqrt(np.sum(q**2)), "orth_pd_vec%d_" % k) # }}} @@ -936,9 +937,51 @@ def tangential_to_xyz(tangential_vec, where=None): for i in range(ambient_dim - 1)) -def project_to_tangential(xyz_vec, which=None): +def project_to_tangential(xyz_vec, where=None): return tangential_to_xyz( - cse(xyz_to_tangential(xyz_vec, which), which)) + cse(xyz_to_tangential(xyz_vec, where), where)) + + +def n_dot(vec, where=None): + nrm = normal(len(vec), where).as_vector() + + return np.dot(nrm, vec) + + +def n_cross(vec, where=None): + nrm = normal(3, where).as_vector() + + from pytools import levi_civita + from pytools.obj_array import make_obj_array + return make_obj_array([ + sum( + levi_civita((i, j, k)) * nrm[j] * vec[k] + for j in range(3) for k in range(3)) + for i in range(3)]) + + +def div_S(kernel, arg, ambient_dim=None, **kwargs): # noqa: N802 + if ambient_dim is None: + ambient_dim = len(arg) + + return sum(dd_axis(iaxis, ambient_dim, S(kernel, arg[iaxis])) + for iaxis in range(ambient_dim)) + + +def curl_S(kernel, arg, **kwargs): # noqa: N802 + from pytools import levi_civita + from pytools.obj_array import make_obj_array + + return make_obj_array([ + sum( + levi_civita((l, m, n)) * dd_axis(m, 3, S(kernel, arg[n], **kwargs)) + for m in range(3) for n in range(3)) + for l in range(3)]) + +# }}} + + +# {{{ vectorial (non-geometric-algebra) differential operators # }}} -- GitLab From 5a1b6eed77a42917f10490024791a69bc31f0b59 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 19:50:27 -0400 Subject: [PATCH 03/11] Modernize PEC C-C MFIE --- pytential/symbolic/pde/maxwell/__init__.py | 34 +++++++++++++--------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index e9301835..35454aff 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -22,10 +22,16 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from six.moves import range, zip import numpy as np # noqa -from pymbolic.primitives import Variable from pytential import sym +from collections import namedtuple +from functools import partial + +tangential_to_xyz = sym.tangential_to_xyz +xyz_to_tangential = sym.xyz_to_tangential +cse = sym.cse + +from pytools.obj_array import join_fields # {{{ Charge-Current MFIE @@ -45,41 +51,41 @@ class PECAugmentedMFIEOperator: def j_operator(self, loc, Jt): Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") return xyz_to_tangential( - (loc*0.5)*Jxyz-nxcurl_S(self.k, 0, Jxyz)) + (loc*0.5)*Jxyz - sym.nxcurl_S(self.kernel, 0, Jxyz, k=self.k)) def j_rhs(self, Hinc_xyz): - return xyz_to_tangential(n_cross(Hinc_xyz)) + return xyz_to_tangential(sym.n_cross(Hinc_xyz)) def rho_operator(self, loc, rho): - return (loc*0.5)*rho+Sp(self.k, rho) + return (loc*0.5)*rho+sym.Sp(self.kernel, rho, k=self.k) def rho_rhs(self, Jt, Einc_xyz): Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") - return n_dot(Einc_xyz) + 1j*self.k*n_dot(S(self.k, Jxyz)) + return (sym.n_dot(Einc_xyz) + + 1j*self.k*sym.n_dot(sym.S(self.kernel, Jxyz, k=self.k))) def scattered_boundary_field(self, Jt, rho, loc): Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") - A = S(self.kernel, Jxyz, k=self.k) - grad_phi = grad_S(k, rho, 3) + A = sym.S(self.kernel, Jxyz, k=self.k) + grad_phi = sym.grad(3, sym.S(self.kernel, rho, k=self.k)) # use - n x n x v = v_tangential - E_scat = 1j*k*A - grad_phi + 0.5*loc*rho - H_scat = curl_S_volume(k, Jxyz) + (loc*0.5)*Jxyz + E_scat = 1j*self.k*A - grad_phi + 0.5*loc*rho + H_scat = sym.curl_S(self.kernel, Jxyz, k=self.k) + (loc*0.5)*Jxyz - from pytools.obj_array import join_fields return join_fields(E_scat, H_scat) def scattered_volume_field(self, Jt, rho): Jxyz = sym.cse(sym.tangential_to_xyz(Jt), "Jxyz") A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=None) - #grad_phi = grad_S(self.kernel, rho, 3, k=self.k) + grad_phi = sym.grad(3, sym.S(self.kernel, rho, 3, k=self.k)) - E_scat = 1j*self.k*A# - grad_phi - #H_scat = curl_S_volume(k, Jxyz) + E_scat = 1j*self.k*A - grad_phi + H_scat = sym.curl_S(self.kernel, Jxyz, k=self.k) from pytools.obj_array import join_fields return E_scat #join_fields(E_scat, H_scat) -- GitLab From 4fc1d146c4f119d71569fcb85311e8b9d941aabb Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 19:51:51 -0400 Subject: [PATCH 04/11] Make Mueller operator actually executable, move generalized Debye out of the way --- pytential/symbolic/pde/maxwell/__init__.py | 465 +++--------------- .../symbolic/pde/maxwell/generalized_debye.py | 391 +++++++++++++++ setup.cfg | 3 +- 3 files changed, 449 insertions(+), 410 deletions(-) create mode 100644 pytential/symbolic/pde/maxwell/generalized_debye.py diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 35454aff..1b0e115a 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -87,34 +87,29 @@ class PECAugmentedMFIEOperator: E_scat = 1j*self.k*A - grad_phi H_scat = sym.curl_S(self.kernel, Jxyz, k=self.k) - from pytools.obj_array import join_fields - return E_scat #join_fields(E_scat, H_scat) + return join_fields(E_scat, H_scat) # }}} # {{{ Charge-Current Mueller MFIE -class MuellerAugmentedMFIEOperator: - """Magnetic Field Integral Equation operator, - under the assumption of no surface charges. - """ - - def __init__(self,omg,mu0,eps0,mu1,eps1): +class MuellerAugmentedMFIEOperator(object): + def __init__(self, omega, mus, epss): from sumpy.kernel import HelmholtzKernel self.kernel = HelmholtzKernel(3) - self.omg = omg - self.mu0 = mu0 - self.eps0 = eps0 - self.mu1=mu1 - self.eps1=eps1 - self.k0 = omg*sqrt(eps0*mu0) - self.k1 = omg*sqrt(eps1*mu1) + self.omega = omega + self.mus = mus + self.epss = epss + self.ks = [ + sym.cse(omega*sym.sqrt(eps*mu), "k%d" % i) + for i, (eps, mu) in enumerate(zip(epss, mus))] def make_unknown(self, name): return sym.make_sym_vector(name, 6) - unk_structure = namedtuple(["jt", "rho_e", "mt", "rho_m"] + unk_structure = namedtuple("MuellerUnknowns", ["jt", "rho_e", "mt", "rho_m"]) + def split_unknown(self, unk): return self.unk_structure( jt=unk[:2], @@ -122,416 +117,70 @@ class MuellerAugmentedMFIEOperator: mt=unk[3:5], rho_m=unk[5]) - def augmul_operator(self, unk): + def operator(self, unk): u = self.split_unknown(unk) - + Jxyz = cse(tangential_to_xyz(u.jt), "Jxyz") Mxyz = cse(tangential_to_xyz(u.mt), "Mxyz") - E0 = cse(1j*self.omg*self.mu0*self.eps0*S(self.k0,0,Jxyz) + - self.mu0*gradScross(self.k0,Mxyz)-grad_S(self.k0,u.rho_e,3),"E0") - H0 = cse(-1j*self.omg*self.mu0*self.eps0*S(self.k0,0,Mxyz) + - self.eps0*gradScross(self.k0,Jxyz)+grad_S(self.k0,u.rho_m,3),"H0") - E1 = cse(1j*self.omg*self.mu1*self.eps1*S(self.k1,0,Jxyz) + - self.mu1*gradScross(self.k1,Mxyz)-grad_S(self.k1,u.rho_e,3),"E1") - H1 = cse(-1j*self.omg*self.mu1*self.eps1*S(self.k1,0,Mxyz) + - self.eps1*gradScross(self.k1,Jxyz)+grad_S(self.k1,u.rho_m,3),"H1") - F1 = cse(xyz_to_tangential(n_cross(H1-H0) - + 0.5*(self.eps0+self.eps1)*Jxyz),"F1") - F2 = cse(n_dot(self.eps1*E1-self.eps0*E0) - + 0.5*(self.eps1+self.eps0)*u.rho_e,"F2") - F3 = cse(xyz_to_tangential(n_cross(E1-E0) - + 0.5*(self.mu0+self.mu1)*Mxyz),"F3") - # sign flip included - F4 = cse(-n_dot(self.mu1*H1-self.mu0*H0) - + 0.5*(self.mu1+self.mu0)*u.rho_m,"F4") - - return sym.join_fields([F1,F2,F3,F4]) - - def augmul_rhs(self, Einc_xyz,Hinc_xyz): - return sym.join_fields([xyz_to_tangential(n_cross(Hinc_xyz)), - n_dot(self.eps1*Einc_xyz),xyz_to_tangential(n_cross(Einc_xyz)), - n_dot(-self.mu1*Hinc_xyz)] - - def scattered_volume_field(self, sol): - u = self.split_unknown(sol) - Jxyz = sym.cse(sym.tangential_to_xyz(u.jt), "Jxyz") - Mxyz = sym.cse(sym.tangential_to_xyz(u.mt), "Mxyz") - - E0 = cse(1j*self.omg*self.mu0*self.eps0*S(self.k0,0,Jxyz) + - self.mu0*gradScross(self.k0,Mxyz)-grad_S(self.k0,u.rho_e,3),"E0") - H0 = cse(-1j*self.omg*self.mu0*self.eps0*S(self.k0,0,Mxyz) + - self.eps0*gradScross(self.k0,Jxyz)+grad_S(self.k0,u.rho_m,3),"H0") - - from pytools.obj_array import join_fields - return join_fields(E0, H0) - -# }}} - - -# {{{ generalized Debye - -def _debye_S0_surf_div(nxnxE): - # Note integration by parts: S0 div_\Gamma (X) -> - \int \grad g0 \cdot - - nxnxE_t = cse(xyz_to_tangential(nxnxE), "nxnxE_t") - - return - sum([ - IntGdSource(0, nxnxE_t[i], - ds_direction=make_tangent(i, 3)) - for i in range(3-1)]) - - -class DebyeOperatorBase(object): - """ - See the `wiki page `_ - for detailed documentation. - - j, m are the (non-physical) Debye currents - r, q are (non-physical) Debye charge densities (from the Debye paper) - - r_tilde, defined by: r = surf_lap S_0^2 r_tilde - q_tilde, defined by: q = surf_lap S_0^2 q_tilde - - :class:`InvertingDebyeOperator` below implements the version based - on r and q. - - :class:`NonInvertingDebyeOperator` below implements the version based - on r_tilde and q_tilde. - - A, Q are the Debye vector potentials - phi, psi are the Debye scalar potentials - """ - - def __init__(self, loc_sign, k, invertible, genus=0, - a_cycle_names=None, - b_cycle_names=None, b_spanning_surface_names=None, - harmonic_vector_field_names=None, h_on_spanning_surface_names=None): - self.loc_sign = loc_sign - - self.k = k - self.invertible = invertible - - self.genus = genus - - # {{{ symbolic names for multiply connected case - - if harmonic_vector_field_names is None: - harmonic_vector_field_names = \ - ["hvf%d" % i for i in range(2*genus)] - - self.harmonic_vector_field_names = harmonic_vector_field_names - self.harmonic_vector_field_symbols = [ - make_vector_field(name, 3) - for name in harmonic_vector_field_names] - - if a_cycle_names is None: - a_cycle_names = ["acyc%d" % i for i in range(genus)] - self.a_cycle_names = a_cycle_names - - if b_cycle_names is None: - b_cycle_names = ["bcyc%d" % i for i in range(genus)] - self.b_cycle_names = b_cycle_names - - if b_spanning_surface_names is None: - b_spanning_surface_names = ["span_surf_%d" % i for i in range(genus)] - self.b_spanning_surface_names = b_spanning_surface_names - - if h_on_spanning_surface_names is None: - h_on_spanning_surface_names = ["h_%s" % name - for name in b_spanning_surface_names] - self.h_on_spanning_surface_names = h_on_spanning_surface_names - - self.h_on_spanning_surface_symbols = [ - make_vector_field(name, 3) - for name in h_on_spanning_surface_names] - - # }}} - - def m(self, j): - return n_cross(j) - - def boundary_field_base(self, r, q, j): - k = self.k - - surf_grad_phi = surf_grad_S(self.k, r) - - dpsi_dn = Sp(self.k, q) - dpsi_dn -= self.loc_sign*1/2*q - - m = cse(self.m(j), "m") - - A = cse(S(k, j), "A") - Q = cse(S(k, m), "Q") - - nxnxE = ( - n_cross( - n_cross(1j * k * A) - - - nxcurl_S(k, self.loc_sign, m)) - - # sign: starts out as -grad phi - # nxnx = - tangential component - + surf_grad_phi) - - ndotH = ( - np.dot(make_normal(3), - # no limit terms: normal part of curl A - curl_S_volume(k, j) - - + 1j*k*Q) - - - dpsi_dn) - - return nxnxE, ndotH - def volume_field_base(self, r, q, j): - k = self.k + omega = self.omega + mu0, mu1 = self.mus + eps0, eps1 = self.epss + k0, k1 = self.ks - grad_phi = grad_S(k, r, 3) - grad_psi = grad_S(k, q, 3) + S = partial(sym.S, self.kernel, qbx_forced_limit="avg") + curl_S = partial(sym.curl_S, self.kernel, qbx_forced_limit="avg") + grad = partial(sym.grad, 3) - m = cse(self.m(j), "m") + E0 = sym.cse(1j*omega*mu0*eps0*S(Jxyz, k=k0) + + mu0*curl_S(Mxyz, k=k0) - grad(S(u.rho_e, k=k0)), "E0") + H0 = sym.cse(-1j*omega*mu0*eps0*S(Mxyz, k=k0) + + eps0*curl_S(Jxyz, k=k0) + grad(S(u.rho_m, k=k0)), "H0") + E1 = sym.cse(1j*omega*mu1*eps1*S(Jxyz, k=k1) + + mu1*curl_S(Mxyz, k=k1) - grad(S(u.rho_e, k=k1)), "E1") + H1 = sym.cse(-1j*omega*mu1*eps1*S(Mxyz, k=k1) + + eps1*curl_S(Jxyz, k=k1) + grad(S(u.rho_m, k=k1)), "H1") - A = cse(S(k, j), "A") - Q = cse(S(k, m), "Q") + F1 = (xyz_to_tangential(sym.n_cross(H1-H0) + 0.5*(eps0+eps1)*Jxyz)) + F2 = (sym.n_dot(eps1*E1-eps0*E0) + 0.5*(eps1+eps0)*u.rho_e) + F3 = (xyz_to_tangential(sym.n_cross(E1-E0) + 0.5*(mu0+mu1)*Mxyz)) - E = 1j*k*A - grad_phi - curl_S_volume(k, m) - H = curl_S_volume(k, j) + 1j*k*Q - grad_psi + # sign flip included + F4 = - sym.n_dot(mu1*H1-mu0*H0) + 0.5*(mu1+mu0)*u.rho_m - from pytools.obj_array import join_fields - return join_fields(E, H) - - def integral_equation(self, *args, **kwargs): - nxnxE, ndotH = self.boundary_field(*args) - nxnxE = cse(nxnxE, "nxnxE") - - fix = kwargs.pop("fix", 0) - if kwargs: - raise TypeError("invalid keyword argument(s)") - - from pytools.obj_array import make_obj_array - eh_op = make_obj_array([ - 2*_debye_S0_surf_div(nxnxE), - -ndotH, - ]) + fix + return join_fields(F1, F2, F3, F4) - k = self.k - j = cse(self.j(*args), "j") - m = cse(self.m(j), "m") - A = cse(S(k, j), "A") + def rhs(self, Einc_xyz, Hinc_xyz): + mu1 = self.mus[1] + eps1 = self.epss[1] - E_minus_grad_phi = 1j*k*A - curl_S_volume(k, m) - - from hellskitchen.fmm import DifferenceKernel - from pytools.obj_array import join_fields return join_fields( - eh_op, - # FIXME: These are inefficient. They compute a full volume field, - # but only actually use the line part of it. - - [ - # Grad phi integrated on a loop does not contribute. - LineIntegral(E_minus_grad_phi, a_cycle_name) - for a_cycle_name in self.a_cycle_names - ], - [ - LineIntegral( - # (E(k) - E(0))/(1jk) - # Grad phi integrated on a loop does not contribute. - (1j*k*A - curl_S_volume(DifferenceKernel(k), m)) - /(1j*k), - b_cycle_name) - for b_cycle_name in self.b_cycle_names - ] - ) - - def prepare_rhs(self, e_inc, h_inc): - normal = make_normal(3) - nxnxE = cse(- project_to_tangential(e_inc), "nxnxE") - - e_rhs = 2*cse(_debye_S0_surf_div(nxnxE)) - h_rhs = cse(-np.dot(normal, h_inc)) - - result = [e_rhs, h_rhs] - - result.extend([ - -LineIntegral(h_inc, a_cycle_name) - for a_cycle_name in self.a_cycle_names - ] + [ - Integral( - np.dot(make_normal(3, ssurf_name), h_on_surf), - ssurf_name) - - for ssurf_name, h_on_surf in zip( - self.b_spanning_surface_names, - self.h_on_spanning_surface_symbols)] - ) - - from pytools.obj_array import make_obj_array - return make_obj_array(result) - - - def harmonic_vector_field_current(self, hvf_coefficients): - return sum(hvf_coeff_i * hvf_i - for hvf_coeff_i, hvf_i in - zip(hvf_coefficients, self.harmonic_vector_field_symbols)) - - def cluster_points(self): - return (-1/4, +1/2) - - - - - -class InvertingDebyeOperatorBase(DebyeOperatorBase): - def r_and_q(self, r, q): - return r, q - - def boundary_field(self, r, q, hvf_coefficients): - j = cse(self.j(r, q, hvf_coefficients), "j") - return self.boundary_field_base(r, q, j) - - def volume_field(self, r, q, hvf_coefficients): - j = cse(self.j(r, q, hvf_coefficients), "j") - return self.volume_field_base(r, q, j) - - def integral_equation(self, r_tilde, q_tilde, hvf_coefficients): - fix = 0 - - if self.invertible: - s_ones = cse(S(0,Ones()), "s_ones") - - def inv_rank_one_coeff(u): - return cse(Mean(u)) - - r_coeff = inv_rank_one_coeff(r_tilde) - q_coeff = inv_rank_one_coeff(q_tilde) - - from pytools.obj_array import join_fields - factors = self.cluster_points() - - fix = join_fields( - factors[0]*s_ones*r_coeff, - factors[1]*Ones()*q_coeff, - ) + xyz_to_tangential(sym.n_cross(Hinc_xyz)), + sym.n_dot(eps1*Einc_xyz), + xyz_to_tangential(sym.n_cross(Einc_xyz)), + sym.n_dot(-mu1*Hinc_xyz)) - return DebyeOperatorBase.integral_equation( - self, r_tilde, q_tilde, hvf_coefficients, fix=fix) - - - - - -class InvertingDebyeOperator(InvertingDebyeOperatorBase): - "Debye operator based on r and q." - - def j(self, r, q, hvf_coefficients): - # We would like to solve - # - # surf_lap alpha = f - # - # Let alpha = S^2 sigma. Then solve - # - # surf_lap S^2 sigma = f - # - # instead. (We need the inner S^2 because that's what we can represent.) - - sigma_solve = Variable("sigma") - surf_lap_op = surface_laplacian_S_squared( - sigma_solve, invertibility_scale=1) - sigma_r = cse(IterativeInverse(surf_lap_op, r, "sigma")) - sigma_q = cse(IterativeInverse(surf_lap_op, q, "sigma")) - - surf_grad_alpha_r = cse(surf_grad_S(0, S(0, sigma_r)), "surf_grad_alpha_r") - surf_grad_alpha_q = cse(surf_grad_S(0, S(0, sigma_q)), "surf_grad_alpha_q") - - return ( - 1j * self.k * ( - surf_grad_alpha_r - n_cross(surf_grad_alpha_q)) - + self.harmonic_vector_field_current(hvf_coefficients)) - - -class InvertingSLapSDebyeOperator(InvertingDebyeOperatorBase): - "Debye operator based on r and q." - - def j(self, r, q, hvf_coefficients): - # We would like to solve - # - # surf_lap alpha = f - # - # Let alpha = S sigma. Then solve - # - # S surf_lap S sigma = S f - # - # instead. (We need the inner S and the outer S because that's - # what we can represent--in this version.) - - sigma_solve = Variable("sigma") - surf_lap_op = S_surface_laplacian_S(sigma_solve, 3, invertibility_scale=1) - sigma_r = IterativeInverse(surf_lap_op, S(0, r), "sigma") - sigma_q = IterativeInverse(surf_lap_op, S(0, q), "sigma") - - surf_grad_alpha_r = cse(surf_grad_S(0, sigma_r), "surf_grad_alpha_r") - surf_grad_alpha_q = cse(surf_grad_S(0, sigma_q), "surf_grad_alpha_q") - - return ( - 1j * self.k * ( - surf_grad_alpha_r - n_cross(surf_grad_alpha_q)) - + self.harmonic_vector_field_current(hvf_coefficients)) - - -class NonInvertingDebyeOperator(DebyeOperatorBase): - "Debye operator based on r_tilde and q_tilde." - - def r_and_q(self, r_tilde, q_tilde): - r = cse(surface_laplacian_S_squared(r_tilde), "r") - q = cse(surface_laplacian_S_squared(q_tilde), "q") - return r, q - - def j(self, r_tilde, q_tilde, hvf_coefficients): - s_r_tilde = cse(S(0, r_tilde), "s_r_tilde") - s_q_tilde = cse(S(0, q_tilde), "s_q_tilde") - assert len(hvf_coefficients) == len(self.harmonic_vector_field_symbols) - surf_grad_s2_r_tilde = cse(surf_grad_S(0, s_r_tilde), "surf_grad_s2_r_tilde") - surf_grad_s2_q_tilde = cse(surf_grad_S(0, s_q_tilde), "surf_grad_s2_q_tilde") - - return (1j * self.k * ( - surf_grad_s2_r_tilde - n_cross(surf_grad_s2_q_tilde)) - + self.harmonic_vector_field_current(hvf_coefficients)) - - def boundary_field(self, r_tilde, q_tilde, hvf_coefficients): - r, q = self.r_and_q(r_tilde, q_tilde) - j = cse(self.j(r_tilde, q_tilde, hvf_coefficients), "j") - return self.boundary_field_base(r, q, j) - - def volume_field(self, r_tilde, q_tilde, hvf_coefficients): - r, q = self.r_and_q(r_tilde, q_tilde) - j = cse(self.j(r_tilde, q_tilde, hvf_coefficients), "j") - return self.volume_field_base(r, q, j) + def representation(self, i, sol): + u = self.split_unknown(sol) + Jxyz = sym.cse(sym.tangential_to_xyz(u.jt), "Jxyz") + Mxyz = sym.cse(sym.tangential_to_xyz(u.mt), "Mxyz") - def integral_equation(self, r_tilde, q_tilde, hvf_coefficients): - fix = 0 + omega = self.omega + mu = self.mus[i] + eps = self.epss[i] + k = self.ks[i] - if self.invertible: - s_ones = cse(S(0, Ones()), "s_ones") + S = partial(sym.S, self.kernel, qbx_forced_limit=None, k=k) + curl_S = partial(sym.curl_S, self.kernel, qbx_forced_limit="avg", k=k) + grad = partial(sym.grad, 3) - def inv_rank_one_coeff(u): - return cse(Mean(cse(S(0, cse(S(0, u)))))) + E0 = 1j*omega*mu*eps*S(Jxyz) + mu*curl_S(Mxyz) - grad(S(u.rho_e)) + H0 = -1j*omega*mu*eps*S(Mxyz) + eps*curl_S(Jxyz) + grad(S(u.rho_m)) - r_coeff = inv_rank_one_coeff(r_tilde) - q_coeff = inv_rank_one_coeff(q_tilde) - - from pytools.obj_array import join_fields - factors = self.cluster_points() - - fix = join_fields( - factors[0]*s_ones*(r_coeff), - factors[1]*Ones()*(q_coeff), - ) - - return DebyeOperatorBase.integral_equation( - self, r_tilde, q_tilde, hvf_coefficients, fix=fix) + from pytools.obj_array import join_fields + return join_fields(E0, H0) # }}} + # vim: foldmethod=marker diff --git a/pytential/symbolic/pde/maxwell/generalized_debye.py b/pytential/symbolic/pde/maxwell/generalized_debye.py new file mode 100644 index 00000000..1b6ea8ee --- /dev/null +++ b/pytential/symbolic/pde/maxwell/generalized_debye.py @@ -0,0 +1,391 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2010-2013 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. +""" + + +# {{{ generalized Debye + +def _debye_S0_surf_div(nxnxE): + # Note integration by parts: S0 div_\Gamma (X) -> - \int \grad g0 \cdot + + nxnxE_t = cse(xyz_to_tangential(nxnxE), "nxnxE_t") + + return - sum([ + IntGdSource(0, nxnxE_t[i], + ds_direction=make_tangent(i, 3)) + for i in range(3-1)]) + + +class DebyeOperatorBase(object): + """ + See the `wiki page `_ + for detailed documentation. + + j, m are the (non-physical) Debye currents + r, q are (non-physical) Debye charge densities (from the Debye paper) + + r_tilde, defined by: r = surf_lap S_0^2 r_tilde + q_tilde, defined by: q = surf_lap S_0^2 q_tilde + + :class:`InvertingDebyeOperator` below implements the version based + on r and q. + + :class:`NonInvertingDebyeOperator` below implements the version based + on r_tilde and q_tilde. + + A, Q are the Debye vector potentials + phi, psi are the Debye scalar potentials + """ + + def __init__(self, loc_sign, k, invertible, genus=0, + a_cycle_names=None, + b_cycle_names=None, b_spanning_surface_names=None, + harmonic_vector_field_names=None, h_on_spanning_surface_names=None): + self.loc_sign = loc_sign + + self.k = k + self.invertible = invertible + + self.genus = genus + + # {{{ symbolic names for multiply connected case + + if harmonic_vector_field_names is None: + harmonic_vector_field_names = \ + ["hvf%d" % i for i in range(2*genus)] + + self.harmonic_vector_field_names = harmonic_vector_field_names + self.harmonic_vector_field_symbols = [ + make_vector_field(name, 3) + for name in harmonic_vector_field_names] + + if a_cycle_names is None: + a_cycle_names = ["acyc%d" % i for i in range(genus)] + self.a_cycle_names = a_cycle_names + + if b_cycle_names is None: + b_cycle_names = ["bcyc%d" % i for i in range(genus)] + self.b_cycle_names = b_cycle_names + + if b_spanning_surface_names is None: + b_spanning_surface_names = ["span_surf_%d" % i for i in range(genus)] + self.b_spanning_surface_names = b_spanning_surface_names + + if h_on_spanning_surface_names is None: + h_on_spanning_surface_names = ["h_%s" % name + for name in b_spanning_surface_names] + self.h_on_spanning_surface_names = h_on_spanning_surface_names + + self.h_on_spanning_surface_symbols = [ + make_vector_field(name, 3) + for name in h_on_spanning_surface_names] + + # }}} + + def m(self, j): + return n_cross(j) + + def boundary_field_base(self, r, q, j): + k = self.k + + surf_grad_phi = surf_grad_S(self.k, r) + + dpsi_dn = Sp(self.k, q) + dpsi_dn -= self.loc_sign*1/2*q + + m = cse(self.m(j), "m") + + A = cse(S(k, j), "A") + Q = cse(S(k, m), "Q") + + nxnxE = ( + n_cross( + n_cross(1j * k * A) + + - nxcurl_S(k, self.loc_sign, m)) + + # sign: starts out as -grad phi + # nxnx = - tangential component + + surf_grad_phi) + + ndotH = ( + np.dot(make_normal(3), + # no limit terms: normal part of curl A + curl_S_volume(k, j) + + + 1j*k*Q) + + - dpsi_dn) + + return nxnxE, ndotH + + def volume_field_base(self, r, q, j): + k = self.k + + grad_phi = grad_S(k, r, 3) + grad_psi = grad_S(k, q, 3) + + m = cse(self.m(j), "m") + + A = cse(S(k, j), "A") + Q = cse(S(k, m), "Q") + + E = 1j*k*A - grad_phi - curl_S_volume(k, m) + H = curl_S_volume(k, j) + 1j*k*Q - grad_psi + + from pytools.obj_array import join_fields + return join_fields(E, H) + + def integral_equation(self, *args, **kwargs): + nxnxE, ndotH = self.boundary_field(*args) + nxnxE = cse(nxnxE, "nxnxE") + + fix = kwargs.pop("fix", 0) + if kwargs: + raise TypeError("invalid keyword argument(s)") + + from pytools.obj_array import make_obj_array + eh_op = make_obj_array([ + 2*_debye_S0_surf_div(nxnxE), + -ndotH, + ]) + fix + + k = self.k + j = cse(self.j(*args), "j") + m = cse(self.m(j), "m") + A = cse(S(k, j), "A") + + E_minus_grad_phi = 1j*k*A - curl_S_volume(k, m) + + from hellskitchen.fmm import DifferenceKernel + from pytools.obj_array import join_fields + return join_fields( + eh_op, + # FIXME: These are inefficient. They compute a full volume field, + # but only actually use the line part of it. + + [ + # Grad phi integrated on a loop does not contribute. + LineIntegral(E_minus_grad_phi, a_cycle_name) + for a_cycle_name in self.a_cycle_names + ], + [ + LineIntegral( + # (E(k) - E(0))/(1jk) + # Grad phi integrated on a loop does not contribute. + (1j*k*A - curl_S_volume(DifferenceKernel(k), m)) + /(1j*k), + b_cycle_name) + for b_cycle_name in self.b_cycle_names + ] + ) + + def prepare_rhs(self, e_inc, h_inc): + normal = make_normal(3) + nxnxE = cse(- project_to_tangential(e_inc), "nxnxE") + + e_rhs = 2*cse(_debye_S0_surf_div(nxnxE)) + h_rhs = cse(-np.dot(normal, h_inc)) + + result = [e_rhs, h_rhs] + + result.extend([ + -LineIntegral(h_inc, a_cycle_name) + for a_cycle_name in self.a_cycle_names + ] + [ + Integral( + np.dot(make_normal(3, ssurf_name), h_on_surf), + ssurf_name) + + for ssurf_name, h_on_surf in zip( + self.b_spanning_surface_names, + self.h_on_spanning_surface_symbols)] + ) + + from pytools.obj_array import make_obj_array + return make_obj_array(result) + + + def harmonic_vector_field_current(self, hvf_coefficients): + return sum(hvf_coeff_i * hvf_i + for hvf_coeff_i, hvf_i in + zip(hvf_coefficients, self.harmonic_vector_field_symbols)) + + def cluster_points(self): + return (-1/4, +1/2) + + + + + +class InvertingDebyeOperatorBase(DebyeOperatorBase): + def r_and_q(self, r, q): + return r, q + + def boundary_field(self, r, q, hvf_coefficients): + j = cse(self.j(r, q, hvf_coefficients), "j") + return self.boundary_field_base(r, q, j) + + def volume_field(self, r, q, hvf_coefficients): + j = cse(self.j(r, q, hvf_coefficients), "j") + return self.volume_field_base(r, q, j) + + def integral_equation(self, r_tilde, q_tilde, hvf_coefficients): + fix = 0 + + if self.invertible: + s_ones = cse(S(0,Ones()), "s_ones") + + def inv_rank_one_coeff(u): + return cse(Mean(u)) + + r_coeff = inv_rank_one_coeff(r_tilde) + q_coeff = inv_rank_one_coeff(q_tilde) + + from pytools.obj_array import join_fields + factors = self.cluster_points() + + fix = join_fields( + factors[0]*s_ones*r_coeff, + factors[1]*Ones()*q_coeff, + ) + + return DebyeOperatorBase.integral_equation( + self, r_tilde, q_tilde, hvf_coefficients, fix=fix) + + + + + +class InvertingDebyeOperator(InvertingDebyeOperatorBase): + "Debye operator based on r and q." + + def j(self, r, q, hvf_coefficients): + # We would like to solve + # + # surf_lap alpha = f + # + # Let alpha = S^2 sigma. Then solve + # + # surf_lap S^2 sigma = f + # + # instead. (We need the inner S^2 because that's what we can represent.) + + sigma_solve = Variable("sigma") + surf_lap_op = surface_laplacian_S_squared( + sigma_solve, invertibility_scale=1) + sigma_r = cse(IterativeInverse(surf_lap_op, r, "sigma")) + sigma_q = cse(IterativeInverse(surf_lap_op, q, "sigma")) + + surf_grad_alpha_r = cse(surf_grad_S(0, S(0, sigma_r)), "surf_grad_alpha_r") + surf_grad_alpha_q = cse(surf_grad_S(0, S(0, sigma_q)), "surf_grad_alpha_q") + + return ( + 1j * self.k * ( + surf_grad_alpha_r - n_cross(surf_grad_alpha_q)) + + self.harmonic_vector_field_current(hvf_coefficients)) + + +class InvertingSLapSDebyeOperator(InvertingDebyeOperatorBase): + "Debye operator based on r and q." + + def j(self, r, q, hvf_coefficients): + # We would like to solve + # + # surf_lap alpha = f + # + # Let alpha = S sigma. Then solve + # + # S surf_lap S sigma = S f + # + # instead. (We need the inner S and the outer S because that's + # what we can represent--in this version.) + + sigma_solve = Variable("sigma") + surf_lap_op = S_surface_laplacian_S(sigma_solve, 3, invertibility_scale=1) + sigma_r = IterativeInverse(surf_lap_op, S(0, r), "sigma") + sigma_q = IterativeInverse(surf_lap_op, S(0, q), "sigma") + + surf_grad_alpha_r = cse(surf_grad_S(0, sigma_r), "surf_grad_alpha_r") + surf_grad_alpha_q = cse(surf_grad_S(0, sigma_q), "surf_grad_alpha_q") + + return ( + 1j * self.k * ( + surf_grad_alpha_r - n_cross(surf_grad_alpha_q)) + + self.harmonic_vector_field_current(hvf_coefficients)) + + +class NonInvertingDebyeOperator(DebyeOperatorBase): + "Debye operator based on r_tilde and q_tilde." + + def r_and_q(self, r_tilde, q_tilde): + r = cse(surface_laplacian_S_squared(r_tilde), "r") + q = cse(surface_laplacian_S_squared(q_tilde), "q") + return r, q + + def j(self, r_tilde, q_tilde, hvf_coefficients): + s_r_tilde = cse(S(0, r_tilde), "s_r_tilde") + s_q_tilde = cse(S(0, q_tilde), "s_q_tilde") + assert len(hvf_coefficients) == len(self.harmonic_vector_field_symbols) + surf_grad_s2_r_tilde = cse(surf_grad_S(0, s_r_tilde), "surf_grad_s2_r_tilde") + surf_grad_s2_q_tilde = cse(surf_grad_S(0, s_q_tilde), "surf_grad_s2_q_tilde") + + return (1j * self.k * ( + surf_grad_s2_r_tilde - n_cross(surf_grad_s2_q_tilde)) + + self.harmonic_vector_field_current(hvf_coefficients)) + + def boundary_field(self, r_tilde, q_tilde, hvf_coefficients): + r, q = self.r_and_q(r_tilde, q_tilde) + j = cse(self.j(r_tilde, q_tilde, hvf_coefficients), "j") + return self.boundary_field_base(r, q, j) + + def volume_field(self, r_tilde, q_tilde, hvf_coefficients): + r, q = self.r_and_q(r_tilde, q_tilde) + j = cse(self.j(r_tilde, q_tilde, hvf_coefficients), "j") + return self.volume_field_base(r, q, j) + + def integral_equation(self, r_tilde, q_tilde, hvf_coefficients): + fix = 0 + + if self.invertible: + s_ones = cse(S(0, Ones()), "s_ones") + + def inv_rank_one_coeff(u): + return cse(Mean(cse(S(0, cse(S(0, u)))))) + + r_coeff = inv_rank_one_coeff(r_tilde) + q_coeff = inv_rank_one_coeff(q_tilde) + + from pytools.obj_array import join_fields + factors = self.cluster_points() + + fix = join_fields( + factors[0]*s_ones*(r_coeff), + factors[1]*Ones()*(q_coeff), + ) + + return DebyeOperatorBase.integral_equation( + self, r_tilde, q_tilde, hvf_coefficients, fix=fix) + +# }}} + diff --git a/setup.cfg b/setup.cfg index 3662361d..ec49074b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,7 +1,6 @@ [flake8] -ignore = E126,E127,E128,E123,E226,E241,E242,E265,E402,W503 +ignore = E126,E127,E128,E123,E226,E241,E242,E265,E402,W503,N803,N806,N802 max-line-length=85 exclude= - pytential/symbolic/pde/maxwell/__init__.py, pytential/symbolic/old_diffop_primitives.py, -- GitLab From c551a027f111008becb264180c2b0b9d03925ad6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 19:54:25 -0400 Subject: [PATCH 05/11] Update Maxwell example to exercise C-C Mueller symbolics --- examples/maxwell.py | 48 ++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 43 insertions(+), 5 deletions(-) diff --git a/examples/maxwell.py b/examples/maxwell.py index b3f016a8..8230aecc 100644 --- a/examples/maxwell.py +++ b/examples/maxwell.py @@ -56,13 +56,51 @@ def main(): qbx = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order, fmm_order=qbx_order + 10, fmm_backend="fmmlib") - from pytential.symbolic.pde.maxwell import PECAugmentedMFIEOperator - pde_op = PECAugmentedMFIEOperator() + from pytential.symbolic.pde.maxwell import MuellerAugmentedMFIEOperator + pde_op = MuellerAugmentedMFIEOperator( + omega=3, + epss=[1.5, 1], + mus=[1, 1], + ) from pytential import bind, sym - rho_sym = sym.var("rho") - # jt_sym = sym.make_sym_vector("jt", 2) - # repr_op = pde_op.scattered_volume_field(jt_sym, rho_sym) + unk = pde_op.make_unknown("sigma") + sym_operator = pde_op.operator(unk) + sym_rhs = pde_op.rhs( + sym.make_sym_vector("Einc", 3), + sym.make_sym_vector("Hinc", 3)) + sym_repr = pde_op.representation(0, unk) + + if 1: + expr = sym_repr + print(sym.pretty(expr)) + + print("#"*80) + bound_op = bind(qbx, expr) + print(bound_op.code) + 1/0 + + if 0: + bvp_rhs = bind(qbx, sym_rhs)(queue, + Einc=make_obj_array([ + ... + ]), + Hinc=make_obj_array([ + ... + ]), + ) + + bound_op = bind(qbx, sym_operator) + + from pytential.solve import gmres + gmres_result = gmres( + bound_op.scipy_op(queue, "sigma", dtype=np.complex128, k=k), + bvp_rhs, tol=1e-8, progress=True, + stall_iterations=0, + hard_failure=True) + + sigma = gmres_result.solution + # {{{ make a density -- GitLab From 1e1d2d36a7482975b82f28af96d5ba3940d72d37 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 20:11:21 -0400 Subject: [PATCH 06/11] Remove moved diff op primitives from diff op dumpster [ci skip] --- pytential/symbolic/old_diffop_primitives.py | 31 --------------------- 1 file changed, 31 deletions(-) diff --git a/pytential/symbolic/old_diffop_primitives.py b/pytential/symbolic/old_diffop_primitives.py index 23028e14..66f85d8d 100644 --- a/pytential/symbolic/old_diffop_primitives.py +++ b/pytential/symbolic/old_diffop_primitives.py @@ -39,21 +39,6 @@ def surf_grad_S(kernel, arg, dim): return project_to_tangential(cse(grad_S(kernel, arg, dim))) -def div_S_volume(kernel, arg): - return sum(IntGdTarget(kernel, arg_n, n) for n, arg_n in enumerate(arg)) - - -def curl_S_volume(kernel, arg): - from pytools import levi_civita - from pytools.obj_array import make_obj_array - - return make_obj_array([ - sum( - levi_civita((l, m, n)) * IntGdTarget(kernel, arg[n], m) - for m in range(3) for n in range(3)) - for l in range(3)]) - - def curl_curl_S_volume(k, arg): # By vector identity, this is grad div S volume + k^2 S_k(arg), # since S_k(arg) satisfies a Helmholtz equation. @@ -147,22 +132,6 @@ def project_to_tangential(xyz_vec, which=None): cse(xyz_to_tangential(xyz_vec, which), which)) -def n_dot(vec, which=None): - return np.dot(normal(len(vec), which), vec) - - -def n_cross(vec, which=None): - nrm = normal(3, which) - - from pytools import levi_civita - from pytools.obj_array import make_obj_array - return make_obj_array([ - sum( - levi_civita((i, j, k)) * nrm[j] * vec[k] - for j in range(3) for k in range(3)) - for i in range(3)]) - - def surf_n_cross(tangential_vec): assert len(tangential_vec) == 2 from pytools.obj_array import make_obj_array -- GitLab From e9e7955a3984b4cc3e3ea2705db34403c8b9edd0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 20:20:51 -0400 Subject: [PATCH 07/11] Fix eval location for the C-C Mueller representation --- pytential/symbolic/pde/maxwell/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 1b0e115a..dd91fab7 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -171,7 +171,7 @@ class MuellerAugmentedMFIEOperator(object): k = self.ks[i] S = partial(sym.S, self.kernel, qbx_forced_limit=None, k=k) - curl_S = partial(sym.curl_S, self.kernel, qbx_forced_limit="avg", k=k) + curl_S = partial(sym.curl_S, self.kernel, qbx_forced_limit=None, k=k) grad = partial(sym.grad, 3) E0 = 1j*omega*mu*eps*S(Jxyz) + mu*curl_S(Mxyz) - grad(S(u.rho_e)) -- GitLab From eb7c5005fa13118ccb4c4c45bfd09fef7e4d2f0b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 22:07:36 -0400 Subject: [PATCH 08/11] Ignore generalized Debye in Flake8 --- setup.cfg | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.cfg b/setup.cfg index ec49074b..ccf2b8a3 100644 --- a/setup.cfg +++ b/setup.cfg @@ -3,4 +3,5 @@ ignore = E126,E127,E128,E123,E226,E241,E242,E265,E402,W503,N803,N806,N802 max-line-length=85 exclude= pytential/symbolic/old_diffop_primitives.py, + pytential/symbolic/pde/maxwell/generalized_debye.py, -- GitLab From 4e66b69f2cf06e00d9ab5e5e275a3957a4f01efc Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 7 Jun 2017 12:12:21 -0400 Subject: [PATCH 09/11] Restore functionality to the potential visualizer in examples/maxwell.py [ci skip] --- examples/maxwell.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/maxwell.py b/examples/maxwell.py index 8230aecc..a76c3459 100644 --- a/examples/maxwell.py +++ b/examples/maxwell.py @@ -140,6 +140,9 @@ def main(): qbx_stick_out = qbx.copy(target_stick_out_factor=0.1) from pytential.target import PointsTarget from pytential.qbx import QBXTargetAssociationFailedException + + rho_sym = sym.var("rho") + try: fld_in_vol = bind( (qbx_stick_out, PointsTarget(fplot.points)), -- GitLab From 9d9123d9a614ee2ddceb9f1bbeaa00d55c9f4666 Mon Sep 17 00:00:00 2001 From: Manas Rachh Date: Wed, 7 Jun 2017 22:27:07 -0400 Subject: [PATCH 10/11] Updating maxwell.py --- examples/maxwell.py | 177 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 142 insertions(+), 35 deletions(-) diff --git a/examples/maxwell.py b/examples/maxwell.py index 8230aecc..9b11268f 100644 --- a/examples/maxwell.py +++ b/examples/maxwell.py @@ -18,7 +18,6 @@ queue = cl.CommandQueue(cl_ctx) target_order = 5 qbx_order = 3 -nelements = 60 mode_nr = 0 h = 1 @@ -27,6 +26,7 @@ k = 4 def main(): + # cl.array.to_device(queue, numpy_array) from meshmode.mesh.io import generate_gmsh, FileSource mesh = generate_gmsh( FileSource("ellipsoid.step"), 2, order=2, @@ -58,9 +58,9 @@ def main(): from pytential.symbolic.pde.maxwell import MuellerAugmentedMFIEOperator pde_op = MuellerAugmentedMFIEOperator( - omega=3, - epss=[1.5, 1], - mus=[1, 1], + omega=0.4, + epss=[1.4, 1.0], + mus=[1.2, 1.0], ) from pytential import bind, sym @@ -76,49 +76,156 @@ def main(): print(sym.pretty(expr)) print("#"*80) - bound_op = bind(qbx, expr) + from pytential.target import PointsTarget + + tgt_points=np.zeros((3,1)) + tgt_points[0,0] = 100 + tgt_points[1,0] = -200 + tgt_points[2,0] = 300 + + bound_op = bind((qbx, PointsTarget(tgt_points)), expr) print(bound_op.code) - 1/0 - if 0: - bvp_rhs = bind(qbx, sym_rhs)(queue, - Einc=make_obj_array([ - ... - ]), - Hinc=make_obj_array([ - ... - ]), - ) + if 1: + + def green3e(x,y,z,source,strength,k): + # electric field corresponding to dyadic green's function + # due to monochromatic electric dipole located at "source". + # "strength" is the the intensity of the dipole. + # E = (I + Hess)(exp(ikr)/r) dot (strength) + # + dx = x - source[0] + dy = y - source[1] + dz = z - source[2] + rr = np.sqrt(dx**2 + dy**2 + dz**2) + + fout = np.exp(1j*k*rr)/rr + evec = fout*strength + qmat = np.zeros((3,3),dtype=np.complex128) + + qmat[0,0]=(2*dx**2-dy**2-dz**2)*(1-1j*k*rr) + qmat[1,1]=(2*dy**2-dz**2-dx**2)*(1-1j*k*rr) + qmat[2,2]=(2*dz**2-dx**2-dy**2)*(1-1j*k*rr) + + qmat[0,0]=qmat[0,0]+(-k**2*dx**2*rr**2) + qmat[1,1]=qmat[1,1]+(-k**2*dy**2*rr**2) + qmat[2,2]=qmat[2,2]+(-k**2*dz**2*rr**2) + + qmat[0,1]=(3-k**2*rr**2-3*1j*k*rr)*(dx*dy) + qmat[1,2]=(3-k**2*rr**2-3*1j*k*rr)*(dy*dz) + qmat[2,0]=(3-k**2*rr**2-3*1j*k*rr)*(dz*dx) + + qmat[1,0]=qmat[0,1] + qmat[2,1]=qmat[1,2] + qmat[0,2]=qmat[2,0] + + fout=np.exp(1j*k*rr)/rr**5/k**2 + + fvec = fout*np.dot(qmat,strength) + evec = evec + fvec + return evec + + def green3m(x,y,z,source,strength,k): + # magnetic field corresponding to dyadic green's function + # due to monochromatic electric dipole located at "source". + # "strength" is the the intensity of the dipole. + # H = curl((I + Hess)(exp(ikr)/r) dot (strength)) = + # strength \cross \grad (exp(ikr)/r) + # + dx = x - source[0] + dy = y - source[1] + dz = z - source[2] + rr = np.sqrt(dx**2 + dy**2 + dz**2) + + fout=(1-1j*k*rr)*np.exp(1j*k*rr)/rr**3 + fvec = np.zeros(3,dtype=np.complex128) + fvec[0] = fout*dx + fvec[1] = fout*dy + fvec[2] = fout*dz + + hvec = np.cross(strength,fvec) + + return hvec + + def dipole3e(x,y,z,source,strength,k): + # + # evalaute electric and magnetic field due + # to monochromatic electric dipole located at "source" + # with intensity "strength" + + evec = green3e(x,y,z,source,strength,k) + evec = evec*1j*k + hvec = green3m(x,y,z,source,strength,k) + return evec,hvec + + def dipole3m(x,y,z,source,strength,k): + # + # evalaute electric and magnetic field due + # to monochromatic magnetic dipole located at "source" + # with intensity "strength" + evec = green3m(x,y,z,source,strength,k) + hvec = green3e(x,y,z,source,strength,k) + hvec = -hvec*1j*k + return evec,hvec + + + def dipole3eall(x,y,z,sources,strengths,k): + ns = len(strengths) + evec = np.zeros(3,dtype=np.complex128) + hvec = np.zeros(3,dtype=np.complex128) + + for i in range(ns): + evect,hvect = dipole3e(x,y,z,sources[i],strengths[i],k) + evec = evec + evect + hvec = hvec + hvect + + nodes = density_discr.nodes().with_queue(queue).get() + source = [0.01,-0.03,0.02] +# source = cl.array.to_device(queue,np.zeros(3)) +# source[0] = 0.01 +# source[1] =-0.03 +# source[2] = 0.02 + strength = np.ones(3) + +# evec = cl.array.to_device(queue,np.zeros((3,len(nodes[0])),dtype=np.complex128)) +# hvec = cl.array.to_device(queue,np.zeros((3,len(nodes[0])),dtype=np.complex128)) + + evec = np.zeros((3,len(nodes[0])),dtype=np.complex128) + hvec = np.zeros((3,len(nodes[0])),dtype=np.complex128) + for i in range(len(nodes[0])): + evec[:,i],hvec[:,i] = dipole3e(nodes[0][i],nodes[1][i],nodes[2][i],source,strength,k) + print(np.shape(hvec)) + print(type(evec)) + print(type(hvec)) + + evec = cl.array.to_device(queue,evec) + hvec = cl.array.to_device(queue,hvec) + + bvp_rhs = bind(qbx, sym_rhs)(queue,Einc=evec,Hinc=hvec) + print(np.shape(bvp_rhs)) + print(type(bvp_rhs)) +# print(bvp_rhs) + 1/-1 bound_op = bind(qbx, sym_operator) from pytential.solve import gmres - gmres_result = gmres( + if 0: + gmres_result = gmres( bound_op.scipy_op(queue, "sigma", dtype=np.complex128, k=k), bvp_rhs, tol=1e-8, progress=True, stall_iterations=0, hard_failure=True) - sigma = gmres_result.solution - + sigma = gmres_result.solution - # {{{ make a density - - nodes = density_discr.nodes().with_queue(queue) - - angle = cl.clmath.atan2(nodes[1], nodes[0]) - - sigma = cl.clmath.cos(mode_nr*angle) - if 0: - sigma = 0*angle - from random import randrange - for i in range(5): - sigma[randrange(len(sigma))] = 1 - - sigma = sigma.astype(np.complex128) - - jt = sym.make_obj_array([sigma, sigma]) - rho = sigma + fld_at_tgt = bind((qbx, PointsTarget(tgt_points)), sym_repr)(queue, + sigma=bvp_rhs,k=k) + fld_at_tgt = np.array([ + fi.get() for fi in fld_at_tgt + ]) + print(fld_at_tgt) + 1/0 # }}} -- GitLab From 524a4bd4645a26a33636cc0d87aaf4d5a0b5717a Mon Sep 17 00:00:00 2001 From: Manas Rachh Date: Wed, 14 Jun 2017 15:00:41 -0400 Subject: [PATCH 11/11] updating maxwell stuff and including sphere tester in test_layer_pot.py [ci skip] --- examples/maxwell_sphere.py | 302 +++++++++++++++++++++ pytential/symbolic/pde/maxwell/__init__.py | 8 +- test/test_layer_pot.py | 37 +++ 3 files changed, 343 insertions(+), 4 deletions(-) create mode 100644 examples/maxwell_sphere.py diff --git a/examples/maxwell_sphere.py b/examples/maxwell_sphere.py new file mode 100644 index 00000000..d0cf8b3d --- /dev/null +++ b/examples/maxwell_sphere.py @@ -0,0 +1,302 @@ +from __future__ import division +import numpy as np +import pyopencl as cl +from sumpy.visualization import FieldPlotter +#from mayavi import mlab +from sumpy.kernel import one_kernel_2d, LaplaceKernel, HelmholtzKernel # noqa + +import faulthandler +from six.moves import range +faulthandler.enable() + +import logging +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + +cl_ctx = cl.create_some_context() +queue = cl.CommandQueue(cl_ctx) + +target_order = 5 +qbx_order = 3 +mode_nr = 0 + +h = 1 + +k = 4 + + +def main(): + # cl.array.to_device(queue, numpy_array) + from meshmode.mesh.io import generate_gmsh, FileSource + from meshmode.mesh.generation import generate_icosphere + from meshmode.mesh.refinement import Refiner + mesh = generate_icosphere(1,target_order) + + refinement_increment = 1 + refiner = Refiner(mesh) + for i in range(refinement_increment): + flags = np.ones(mesh.nelements, dtype=bool) + refiner.refine(flags) + mesh = refiner.get_current_mesh() + + + from meshmode.mesh.processing import perform_flips + # Flip elements--gmsh generates inside-out geometry. + mesh = perform_flips(mesh, np.ones(mesh.nelements)) + + print("%d elements" % mesh.nelements) + + from meshmode.mesh.processing import find_bounding_box + bbox_min, bbox_max = find_bounding_box(mesh) + bbox_center = 0.5*(bbox_min+bbox_max) + bbox_size = max(bbox_max-bbox_min) / 2 + + logger.info("%d elements" % mesh.nelements) + + from pytential.qbx import QBXLayerPotentialSource + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + density_discr = Discretization( + cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + qbx = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order, + fmm_order=False, fmm_backend="fmmlib") + + from pytential.symbolic.pde.maxwell import MuellerAugmentedMFIEOperator + pde_op = MuellerAugmentedMFIEOperator( + omega=1.0, + epss=[1.0, 1.0], + mus=[1.0, 1.0], + ) + from pytential import bind, sym + + unk = pde_op.make_unknown("sigma") + sym_operator = pde_op.operator(unk) + sym_rhs = pde_op.rhs( + sym.make_sym_vector("Einc", 3), + sym.make_sym_vector("Hinc", 3)) + sym_repr = pde_op.representation(1, unk) + + if 1: + expr = sym_repr + print(sym.pretty(expr)) + + print("#"*80) + from pytential.target import PointsTarget + + tgt_points=np.zeros((3,1)) + tgt_points[0,0] = 100 + tgt_points[1,0] = -200 + tgt_points[2,0] = 300 + + bound_op = bind((qbx, PointsTarget(tgt_points)), expr) + print(bound_op.code) + + if 1: + + def green3e(x,y,z,source,strength,k): + # electric field corresponding to dyadic green's function + # due to monochromatic electric dipole located at "source". + # "strength" is the the intensity of the dipole. + # E = (I + Hess)(exp(ikr)/r) dot (strength) + # + dx = x - source[0] + dy = y - source[1] + dz = z - source[2] + rr = np.sqrt(dx**2 + dy**2 + dz**2) + + fout = np.exp(1j*k*rr)/rr + evec = fout*strength + qmat = np.zeros((3,3),dtype=np.complex128) + + qmat[0,0]=(2*dx**2-dy**2-dz**2)*(1-1j*k*rr) + qmat[1,1]=(2*dy**2-dz**2-dx**2)*(1-1j*k*rr) + qmat[2,2]=(2*dz**2-dx**2-dy**2)*(1-1j*k*rr) + + qmat[0,0]=qmat[0,0]+(-k**2*dx**2*rr**2) + qmat[1,1]=qmat[1,1]+(-k**2*dy**2*rr**2) + qmat[2,2]=qmat[2,2]+(-k**2*dz**2*rr**2) + + qmat[0,1]=(3-k**2*rr**2-3*1j*k*rr)*(dx*dy) + qmat[1,2]=(3-k**2*rr**2-3*1j*k*rr)*(dy*dz) + qmat[2,0]=(3-k**2*rr**2-3*1j*k*rr)*(dz*dx) + + qmat[1,0]=qmat[0,1] + qmat[2,1]=qmat[1,2] + qmat[0,2]=qmat[2,0] + + fout=np.exp(1j*k*rr)/rr**5/k**2 + + fvec = fout*np.dot(qmat,strength) + evec = evec + fvec + return evec + + def green3m(x,y,z,source,strength,k): + # magnetic field corresponding to dyadic green's function + # due to monochromatic electric dipole located at "source". + # "strength" is the the intensity of the dipole. + # H = curl((I + Hess)(exp(ikr)/r) dot (strength)) = + # strength \cross \grad (exp(ikr)/r) + # + dx = x - source[0] + dy = y - source[1] + dz = z - source[2] + rr = np.sqrt(dx**2 + dy**2 + dz**2) + + fout=(1-1j*k*rr)*np.exp(1j*k*rr)/rr**3 + fvec = np.zeros(3,dtype=np.complex128) + fvec[0] = fout*dx + fvec[1] = fout*dy + fvec[2] = fout*dz + + hvec = np.cross(strength,fvec) + + return hvec + + def dipole3e(x,y,z,source,strength,k): + # + # evalaute electric and magnetic field due + # to monochromatic electric dipole located at "source" + # with intensity "strength" + + evec = green3e(x,y,z,source,strength,k) + evec = evec*1j*k + hvec = green3m(x,y,z,source,strength,k) +# print(hvec) +# print(strength) + return evec,hvec + + def dipole3m(x,y,z,source,strength,k): + # + # evalaute electric and magnetic field due + # to monochromatic magnetic dipole located at "source" + # with intensity "strength" + evec = green3m(x,y,z,source,strength,k) + hvec = green3e(x,y,z,source,strength,k) + hvec = -hvec*1j*k + return evec,hvec + + + def dipole3eall(x,y,z,sources,strengths,k): + ns = len(strengths) + evec = np.zeros(3,dtype=np.complex128) + hvec = np.zeros(3,dtype=np.complex128) + + for i in range(ns): + evect,hvect = dipole3e(x,y,z,sources[i],strengths[i],k) + evec = evec + evect + hvec = hvec + hvect + + nodes = density_discr.nodes().with_queue(queue).get() + source = [0.01,-0.03,0.02] +# source = cl.array.to_device(queue,np.zeros(3)) +# source[0] = 0.01 +# source[1] =-0.03 +# source[2] = 0.02 + strength = np.ones(3) + +# evec = cl.array.to_device(queue,np.zeros((3,len(nodes[0])),dtype=np.complex128)) +# hvec = cl.array.to_device(queue,np.zeros((3,len(nodes[0])),dtype=np.complex128)) + + evec = np.zeros((3,len(nodes[0])),dtype=np.complex128) + hvec = np.zeros((3,len(nodes[0])),dtype=np.complex128) + for i in range(len(nodes[0])): + evec[:,i],hvec[:,i] = dipole3e(nodes[0][i],nodes[1][i],nodes[2][i],source,strength,k) + print(np.shape(hvec)) + print(type(evec)) + print(type(hvec)) + + evec = cl.array.to_device(queue,evec) + hvec = cl.array.to_device(queue,hvec) + + bvp_rhs = bind(qbx, sym_rhs)(queue,Einc=evec,Hinc=hvec) + print(np.shape(bvp_rhs)) + print(type(bvp_rhs)) +# print(bvp_rhs) + 1/-1 + + bound_op = bind(qbx, sym_operator) + + from pytential.solve import gmres + if 1: + gmres_result = gmres( + bound_op.scipy_op(queue, "sigma", dtype=np.complex128, k=k), + bvp_rhs, tol=1e-8, progress=True, + stall_iterations=0, + hard_failure=True) + + sigma = gmres_result.solution + + fld_at_tgt = bind((qbx, PointsTarget(tgt_points)), sym_repr)(queue, + sigma=sigma,k=k) + fld_at_tgt = np.array([ + fi.get() for fi in fld_at_tgt + ]) + print(fld_at_tgt) + fld_exact_e,fld_exact_h = dipole3e(tgt_points[0,0],tgt_points[1,0],tgt_points[2,0],source,strength,k) + print(fld_exact_e) + print(fld_exact_h) + 1/0 + + # }}} + + #mlab.figure(bgcolor=(1, 1, 1)) + if 1: + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, density_discr, target_order) + + bdry_normals = bind(density_discr, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("source.vtu", [ + ("sigma", sigma), + ("bdry_normals", bdry_normals), + ]) + + fplot = FieldPlotter(bbox_center, extent=2*bbox_size, npoints=(150, 150, 1)) + + qbx_stick_out = qbx.copy(target_stick_out_factor=0.1) + from pytential.target import PointsTarget + from pytential.qbx import QBXTargetAssociationFailedException + + rho_sym = sym.var("rho") + + try: + fld_in_vol = bind( + (qbx_stick_out, PointsTarget(fplot.points)), + sym.make_obj_array([ + sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), + qbx_forced_limit=None), + sym.d_dx(3, sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), + qbx_forced_limit=None)), + sym.d_dy(3, sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), + qbx_forced_limit=None)), + sym.d_dz(3, sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), + qbx_forced_limit=None)), + ]) + )(queue, jt=jt, rho=rho, k=k) + except QBXTargetAssociationFailedException as e: + fplot.write_vtk_file( + "failed-targets.vts", + [ + ("failed_targets", e.failed_target_flags.get(queue)) + ]) + raise + + fld_in_vol = sym.make_obj_array( + [fiv.get() for fiv in fld_in_vol]) + + #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) + fplot.write_vtk_file( + "potential.vts", + [ + ("potential", fld_in_vol[0]), + ("grad", fld_in_vol[1:]), + ] + ) + + +if __name__ == "__main__": + main() diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index dd91fab7..717c7253 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -102,7 +102,7 @@ class MuellerAugmentedMFIEOperator(object): self.mus = mus self.epss = epss self.ks = [ - sym.cse(omega*sym.sqrt(eps*mu), "k%d" % i) + sym.cse(omega*(eps*mu)**0.5, "k%d" % i) for i, (eps, mu) in enumerate(zip(epss, mus))] def make_unknown(self, name): @@ -146,7 +146,7 @@ class MuellerAugmentedMFIEOperator(object): F3 = (xyz_to_tangential(sym.n_cross(E1-E0) + 0.5*(mu0+mu1)*Mxyz)) # sign flip included - F4 = - sym.n_dot(mu1*H1-mu0*H0) + 0.5*(mu1+mu0)*u.rho_m + F4 = -sym.n_dot(mu1*H1-mu0*H0) + 0.5*(mu1+mu0)*u.rho_m return join_fields(F1, F2, F3, F4) @@ -174,8 +174,8 @@ class MuellerAugmentedMFIEOperator(object): curl_S = partial(sym.curl_S, self.kernel, qbx_forced_limit=None, k=k) grad = partial(sym.grad, 3) - E0 = 1j*omega*mu*eps*S(Jxyz) + mu*curl_S(Mxyz) - grad(S(u.rho_e)) - H0 = -1j*omega*mu*eps*S(Mxyz) + eps*curl_S(Jxyz) + grad(S(u.rho_m)) + E0 = 1j*k*eps*S(Jxyz) + mu*curl_S(Mxyz) - grad(S(u.rho_e)) + H0 = -1j*k*mu*S(Mxyz) + eps*curl_S(Jxyz) + grad(S(u.rho_m)) from pytools.obj_array import join_fields return join_fields(E0, H0) diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 7e86a96e..947c1464 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -787,6 +787,43 @@ class EllipsoidIntEqTestCase(IntEqTestCase): # kidding? gmres_tol = 1e-5 +class SphereIntEqTestCase(IntEqTestCase): + resolutions = [2, 1] + name = "sphere" + + def get_mesh(self, resolution, target_order): + from meshmode.mesh.io import generate_gmsh, FileSource + from meshmode.mesh.generation import generate_icosphere + from meshmode.mesh.refinement import Refiner + mesh = generate_icosphere(1,target_order) + + refinement_increment = 1 + refiner = Refiner(mesh) + for i in range(refinement_increment): + flags = np.ones(mesh.nelements, dtype=bool) + refiner.refine(flags) + mesh = refiner.get_current_mesh() + + + from meshmode.mesh.processing import perform_flips + # Flip elements--gmsh generates inside-out geometry. + return perform_flips(mesh, np.ones(mesh.nelements)) + + + fmm_backend = "fmmlib" + use_refinement = False + neumann_alpha = 0 # no double layers in FMMlib backend yet + + inner_radius = 0.4 + outer_radius = 5 + + qbx_order = 2 + target_order = qbx_order + check_tangential_deriv = False + + # We're only expecting three digits based on FMM settings. Who are we + # kidding? + gmres_tol = 1e-5 @pytest.mark.parametrize("case", [ EllipseIntEqTestCase(helmholtz_k=helmholtz_k, bc_type=bc_type, -- GitLab