From 6a46fcd19b920bf93d98c2c60b7fd383300aeb9a Mon Sep 17 00:00:00 2001 From: Manas Rachh Date: Tue, 6 Jun 2017 16:18:10 -0400 Subject: [PATCH 01/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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/53] 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 From 24da4d892b4f1f89d44436f2299059aa93161ca5 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 01:03:20 -0500 Subject: [PATCH 12/53] Add MFIE notes --- contrib/notes/.gitignore | 1 + contrib/notes/mfie.tm | 112 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 113 insertions(+) create mode 100644 contrib/notes/.gitignore create mode 100644 contrib/notes/mfie.tm diff --git a/contrib/notes/.gitignore b/contrib/notes/.gitignore new file mode 100644 index 00000000..55acb67b --- /dev/null +++ b/contrib/notes/.gitignore @@ -0,0 +1 @@ +fontconfig diff --git a/contrib/notes/mfie.tm b/contrib/notes/mfie.tm new file mode 100644 index 00000000..ff3d86b7 --- /dev/null +++ b/contrib/notes/mfie.tm @@ -0,0 +1,112 @@ + + + + +<\body> + + + + + We'll do the derivation for the exterior MFIE from the boundary condition + + <\equation*> + n\-H|)>=J. + + + Next, assume =0> because of PEC, thus + + <\equation*> + n\H=J. + + + Then split into incoming and scattered, i.e. + + <\equation*> + n\+H|)>=J. + + + Use the representation =\\A=\\SJ>, + leading to + + <\equation*> + n\H+\SJ|)>|)>=J. + + + Use \Sv|)>|]>=\Sv|)>|]>+loc>: + (where 1> depending on surface side) + + <\equation*> + n\H+\SJ|)>|)>+|2>=J. + + + Rearrange: + + <\equation*> + n\H=|2>-\SJ|)>|)>. + + + The interior MFIE is derived similarly as: + + <\equation*> + n\H=-|2>-\SJ|)>|)>. + + + > postprocessor derivation> + + We'll start with the boundary condition + + <\equation*> + n\-E|)>=\. + + + Again, =0> because of PEC, i.e. + + <\equation*> + n\+E|)>=\. + + + Now use the representation =i*k*A-\\=i*k*SJ-\S\> + and obtain + + <\equation*> + n\E+J-\S\|)>|)>=\. + + + Carrying out the limit, we obtain: + + <\equation*> + n\E+n\J|)>-S\+\=\. + + + Rearrange: + + <\equation*> + n\E+n\J|)>=\+S\. + + + +<\references> + <\collection> + > + > + > + + + +<\auxiliary> + <\collection> + <\associate|toc> + |math-font-series||Augmented + MFIE derivation> |.>>>>|> + + + |MFIE derivation + |.>>>>|> + > + + ||\> + postprocessor derivation |.>>>>|> + > + + + \ No newline at end of file -- GitLab From 6a1fec06fec5ef5b4504ea5c33e92e9454d721a4 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 01:05:50 -0500 Subject: [PATCH 13/53] Minor Maxwell style/doc tweaks --- pytential/symbolic/pde/maxwell/__init__.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 717c7253..0d0f5af3 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -40,7 +40,7 @@ class PECAugmentedMFIEOperator: """Magnetic Field Integral Equation operator, under the assumption of no surface charges. - see notes/mfie.tm + see :file:`contrib/notes/mfie.tm` """ def __init__(self, k=sym.var("k")): @@ -95,6 +95,10 @@ class PECAugmentedMFIEOperator: # {{{ Charge-Current Mueller MFIE class MuellerAugmentedMFIEOperator(object): + """ + ... warning:: currently untested + """ + def __init__(self, omega, mus, epss): from sumpy.kernel import HelmholtzKernel self.kernel = HelmholtzKernel(3) @@ -165,7 +169,7 @@ class MuellerAugmentedMFIEOperator(object): Jxyz = sym.cse(sym.tangential_to_xyz(u.jt), "Jxyz") Mxyz = sym.cse(sym.tangential_to_xyz(u.mt), "Mxyz") - omega = self.omega + # omega = self.omega mu = self.mus[i] eps = self.epss[i] k = self.ks[i] -- GitLab From c1d5e865eb5efb28d46202d891831038efbfb3f9 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 18:00:29 -0500 Subject: [PATCH 14/53] Doc typo fix --- pytential/source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/source.py b/pytential/source.py index f6cabf20..88dd6963 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -54,7 +54,7 @@ class PotentialSource(object): class PointPotentialSource(PotentialSource): """ - ... attributes:: points + ... attribute:: points An :class:`pyopencl.array.Array` of shape ``[ambient_dim, npoints]``. """ -- GitLab From 1197a44cf522396f95a4ceb4e8d5f4c159d8642d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 18:00:46 -0500 Subject: [PATCH 15/53] Add PointPotentialSource.nnodes --- pytential/source.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/pytential/source.py b/pytential/source.py index 88dd6963..7e4f81fe 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -57,6 +57,8 @@ class PointPotentialSource(PotentialSource): ... attribute:: points An :class:`pyopencl.array.Array` of shape ``[ambient_dim, npoints]``. + + .. attribute:: nnodes """ def __init__(self, cl_context, points): @@ -67,6 +69,10 @@ class PointPotentialSource(PotentialSource): def real_dtype(self): return self.points.dtype + @property + def nnodes(self): + return self.points.shape[-1] + @property def complex_dtype(self): return { -- GitLab From 2c0974af2f61de9cb958da86e899cef568ce07b7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 18:15:09 -0500 Subject: [PATCH 16/53] Fix refiner usage/remove flip in SphereIntEqTestCase --- test/test_scalar_int_eq.py | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py index 706dd855..690541a5 100644 --- a/test/test_scalar_int_eq.py +++ b/test/test_scalar_int_eq.py @@ -611,19 +611,12 @@ class SphereIntEqTestCase(IntEqTestCase): def get_mesh(self, resolution, target_order): from meshmode.mesh.generation import generate_icosphere - from meshmode.mesh.refinement import Refiner - mesh = generate_icosphere(1, target_order) + from meshmode.mesh.refinement import refine_uniformly + mesh = refine_uniformly( + generate_icosphere(1, target_order), + resolution) - 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)) + return mesh fmm_backend = "fmmlib" use_refinement = False -- GitLab From 2b02fd96fb1389aed7a033932389f57dc0d0291f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 18:16:39 -0500 Subject: [PATCH 17/53] Expose join_fields in sym --- pytential/symbolic/primitives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 1343919d..0f93210a 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -35,7 +35,7 @@ from pymbolic.geometric_algebra import MultiVector, componentwise from pymbolic.geometric_algebra.primitives import ( # noqa: F401 NablaComponent, DerivativeSource, Derivative as DerivativeBase) from pymbolic.primitives import make_sym_vector # noqa: F401 -from pytools.obj_array import make_obj_array # noqa: F401 +from pytools.obj_array import make_obj_array, join_fields # noqa: F401 from functools import partial -- GitLab From 94e2da43f9fb29cf7105c7bef6531b279f48bbe4 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 18:16:43 -0500 Subject: [PATCH 18/53] Implement generic 'curl', use in curl_S --- pytential/symbolic/primitives.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 0f93210a..12c96d8c 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -1034,20 +1034,19 @@ def div_S(kernel, arg, ambient_dim=None, **kwargs): # noqa: N802 for iaxis in range(ambient_dim)) -def curl_S(kernel, arg, **kwargs): # noqa: N802 +def curl(vec): 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)) + levi_civita((l, m, n)) * dd_axis(m, 3, vec[n]) for m in range(3) for n in range(3)) for l in range(3)]) -# }}} - -# {{{ vectorial (non-geometric-algebra) differential operators +def curl_S(kernel, arg, **kwargs): # noqa: N802 + return curl(S(kernel, arg, **kwargs)) # }}} -- GitLab From 9104438540ef233c2b5de7cbd58471bd2e932b73 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 18:17:22 -0500 Subject: [PATCH 19/53] Use sym.join_fields in Maxwell operators (instead of importing from pytools) --- pytential/symbolic/pde/maxwell/__init__.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 0d0f5af3..93dfffc6 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -31,8 +31,6 @@ 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 @@ -76,18 +74,18 @@ class PECAugmentedMFIEOperator: 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 - return join_fields(E_scat, H_scat) + return sym.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 = sym.grad(3, sym.S(self.kernel, rho, 3, k=self.k)) + grad_phi = sym.grad(3, sym.S(self.kernel, rho, k=self.k)) E_scat = 1j*self.k*A - grad_phi H_scat = sym.curl_S(self.kernel, Jxyz, k=self.k) - return join_fields(E_scat, H_scat) + return sym.join_fields(E_scat, H_scat) # }}} @@ -152,13 +150,13 @@ class MuellerAugmentedMFIEOperator(object): # sign flip included F4 = -sym.n_dot(mu1*H1-mu0*H0) + 0.5*(mu1+mu0)*u.rho_m - return join_fields(F1, F2, F3, F4) + return sym.join_fields(F1, F2, F3, F4) def rhs(self, Einc_xyz, Hinc_xyz): mu1 = self.mus[1] eps1 = self.epss[1] - return join_fields( + return sym.join_fields( xyz_to_tangential(sym.n_cross(Hinc_xyz)), sym.n_dot(eps1*Einc_xyz), xyz_to_tangential(sym.n_cross(Einc_xyz)), @@ -181,8 +179,7 @@ class MuellerAugmentedMFIEOperator(object): 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) + return sym.join_fields(E0, H0) # }}} -- GitLab From b32e0c56e8944cd1deb095b0ed03a17c0d4bfad5 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 18:18:26 -0500 Subject: [PATCH 20/53] Add fledgling MFIE test that has a valid Maxwell solution --- test/test_maxwell.py | 289 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 289 insertions(+) create mode 100644 test/test_maxwell.py diff --git a/test/test_maxwell.py b/test/test_maxwell.py new file mode 100644 index 00000000..f652ae65 --- /dev/null +++ b/test/test_maxwell.py @@ -0,0 +1,289 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2017 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 # noqa +import pyopencl as cl +import pyopencl.clmath # noqa +import pyopencl.clrandom # noqa +import pytest + +from pytential import bind, sym + +from sumpy.visualization import FieldPlotter # noqa +from sumpy.point_calculus import CalculusPatch, frequency_domain_maxwell +from pytential.target import PointsTarget + +import logging +logger = logging.getLogger(__name__) + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + + +class TestCase: + def __init__(self, k, is_interior, resolutions, qbx_order, + fmm_order): + self.k = k + self.is_interior = is_interior + self.resolutions = resolutions + self.qbx_order = qbx_order + self.fmm_order = fmm_order + + +class SphereTestCase(TestCase): + def get_mesh(self, resolution, target_order): + from meshmode.mesh.generation import generate_icosphere + from meshmode.mesh.refinement import refine_uniformly + return refine_uniformly( + generate_icosphere(2, target_order), + resolution) + + def get_test_mesh(self, target_order): + from meshmode.mesh.generation import generate_icosphere + + if self.is_interior: + + return generate_icosphere(5, target_order) + else: + return generate_icosphere(0.5, target_order) + + def get_source(self, queue): + if self.is_interior: + source_ctr = np.array([[0.35, 0.1, 0.15]]).T + else: + source_ctr = np.array([[5, 0.1, 0.15]]).T + + sources = source_ctr + np.random.rand(3, 10) + from pytential.source import PointPotentialSource + return PointPotentialSource( + queue.context, + cl.array.to_device(queue, sources)) + + +tc_int = SphereTestCase(k=10.2, is_interior=True, resolutions=[0, 1], + qbx_order=3, fmm_order=10) +tc_ext = SphereTestCase(k=10.2, is_interior=False, resolutions=[0, 1], + qbx_order=3, fmm_order=10) + + +def get_sym_maxwell_source(kernel, jxyz, k): + # This ensures div A = 0, which is simply a consequence of div curl S=0. + # This means we use the Coulomb gauge to generate this field. + + A = sym.curl(sym.S(kernel, jxyz, k=k, qbx_forced_limit=None)) + + # https://en.wikipedia.org/w/index.php?title=Maxwell%27s_equations&oldid=798940325#Alternative_formulations + return sym.join_fields( + - 1j*k*A, + sym.curl(A)) + + +@pytest.mark.parametrize("case", [ + tc_int, + tc_ext + ]) +def test_mfie_from_source(ctx_getter, case, visualize=False): + logging.basicConfig(level=logging.INFO) + + cl_ctx = ctx_getter() + queue = cl.CommandQueue(cl_ctx) + + pytest.importorskip("pyfmmlib") + pytest.importorskip("scipy") + + np.random.seed(12) + + # {{{ come up with a solution to Maxwell's equations + + j_sym = sym.make_sym_vector("j", 3) + + from pytential.symbolic.pde.maxwell import PECAugmentedMFIEOperator + mfie = PECAugmentedMFIEOperator() + + test_source = case.get_source(queue) + + calc_patch = CalculusPatch(np.array([-2, 0, 0]), h=0.01) + calc_patch_tgt = PointsTarget(calc_patch.points) + + rng = cl.clrandom.PhiloxGenerator(cl_ctx, seed=12) + src_j = rng.normal(queue, (3, test_source.nnodes), dtype=np.float64) + + pde_test_field = bind( + (test_source, calc_patch_tgt), + get_sym_maxwell_source(mfie.kernel, j_sym, mfie.k) + )(queue, j=src_j, k=case.k) + + from sumpy.tools import vector_from_device + pde_test_e = vector_from_device(queue, pde_test_field[:3]) + pde_test_h = vector_from_device(queue, pde_test_field[3:]) + + print([ + calc_patch.norm(x, np.inf) / calc_patch.norm(pde_test_e, np.inf) + for x in frequency_domain_maxwell( + calc_patch, pde_test_e, pde_test_h, case.k)]) + + + 1/0 + + + + + + mfie_ext_op = xyz_to_tangential(mfie.boundary_field(+1, Jxyz)) + mfie_int_op = xyz_to_tangential(mfie.boundary_field(-1, Jxyz)) + xyz_mfie_bdry_vol_op = mfie.boundary_field(0, Jxyz) + xyz_mfie_vol_op = mfie.volume_field(Jxyz) + + n_cross_op = n_cross(make_vector_field("vec", 3)) + + from pytools.convergence import EOCRecorder + + eoc_rec_nxH = EOCRecorder() + eoc_rec_E = EOCRecorder() + eoc_rec_H = EOCRecorder() + + + for resolution in case.resolutions: + scat_mesh = generate_icosphere(2, i) + if is_interior: + tgt_mesh = generate_icosphere(5, 0) + else: + tgt_mesh = generate_icosphere(0.5, 1) + + scat_discr = FlatConstantDiscretization(scat_mesh) + tgt_discr = FlatConstantDiscretization(tgt_mesh) + + def sbind(op): return bind(op, scat_discr, iprec=iprec) + + source_func = gen_field_multipole + + E_scat, H_scat = source_func(k, is_interior, scat_mesh.centroids) + E_tgt_true, H_tgt_true = source_func(k, is_interior, tgt_mesh.centroids) + + nxHinc_xyz_mid = sbind(n_cross_op)(vec=H_scat) + nxH_tgt_true = bind(n_cross_op, tgt_discr, iprec=iprec)( + vec=H_tgt_true) + + # {{{ system solve + + nxHinc_t_mid = sbind( + xyz_to_tangential(make_vector_field("vec", 3)) + )(vec=nxHinc_xyz_mid) + + if is_interior: + mfie_solve_op = mfie_int_op + else: + mfie_solve_op = mfie_ext_op + + Jt_mid = solve_lin_op( + sbind(mfie_solve_op).scipy_op("Jt"), + nxHinc_t_mid) + + # }}} + + # {{{ error in nxH + + nxH_tgt = bind(xyz_mfie_bdry_vol_op, (scat_discr, tgt_discr), iprec=iprec) \ + (Jt=Jt_mid) + rel_err_nxH = (tgt_discr.norm(nxH_tgt_true+nxH_tgt) + / tgt_discr.norm(nxH_tgt_true)) + + # }}} + + # {{{ error in E, H + + EH_tgt = bind(xyz_mfie_vol_op, (scat_discr, tgt_discr), iprec=iprec) \ + (Jt=Jt_mid) + E_tgt = EH_tgt[:3] + H_tgt = EH_tgt[3:] + rel_err_H = (tgt_discr.norm(H_tgt_true-H_tgt) + / tgt_discr.norm(H_tgt_true)) + rel_err_E = (tgt_discr.norm(E_tgt_true-E_tgt) + / tgt_discr.norm(E_tgt_true)) + + # }}} + + print("ERR", h, rel_err_nxH, rel_err_H, rel_err_E) + + eoc_rec_nxH.add_data_point(h, rel_err_nxH) + eoc_rec_H.add_data_point(h, rel_err_H) + eoc_rec_E.add_data_point(h, rel_err_E) + + # {{{ visualization + if write_vis: + from pyvisfile.silo import SiloFile + from hellskitchen.visualization import add_to_silo_file + from pytools.obj_array import oarray_real_copy + + Jxyz_mid = bind( + tangential_to_xyz(make_vector_field("vec", 2)))(vec=Jt_mid) + + silo = SiloFile("mfie-int%s-%d.silo" % (is_interior, i)) + add_to_silo_file(silo, scat_geo, cell_data=[ + ("nxHinc_mid", oarray_real_copy(nxHinc_xyz_mid)), + ("J_mid", oarray_real_copy(Jxyz_mid)), + ]) + add_to_silo_file(silo, tgt_geo, cell_data=[ + ("nxHinc_tgt", oarray_real_copy(nxH_tgt)), + ("nxHinc_tgt_true", oarray_real_copy(nxH_tgt_true)), + ("Hinc_tgt", oarray_real_copy(H_tgt)), + ("Hinc_tgt_true", oarray_real_copy(H_tgt_true)), + ("Einc_tgt", oarray_real_copy(E_tgt)), + ("Einc_tgt_true", oarray_real_copy(E_tgt_true)), + ], mesh_name="tgt_mesh") + silo.close() + + # }}} + + print("--------------------------------------------------------") + print("is_interior=%s" % case.is_interior) + print("--------------------------------------------------------") + for which_eoc, eoc_rec in [ + ("nxH", eoc_rec_nxH), + ("H", eoc_rec_H), + ("E", eoc_rec_E)]: + print(which_eoc) + print(eoc_rec.pretty_print()) + + if len(eoc_rec.history) > 1: + assert eoc_rec.order_estimate() > 0.8 + + + + + + +# You can test individual routines by typing +# $ python test_maxwell.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from py.test.cmdline import main + main([__file__]) + +# vim: fdm=marker -- GitLab From c6a7793f04671bec0c5edaf51f22115eee7a5417 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 23:02:03 -0500 Subject: [PATCH 21/53] Remove 'data prep' debug print from fmmlib interface --- pytential/qbx/fmmlib.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index f146f2b6..ebadacc8 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -352,8 +352,6 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): source_level_start_ibox, source_mpoles_view = \ self.multipole_expansions_view(multipole_exps, isrc_level) - print("par data prep lev %d" % isrc_level) - tgt_icenter_vec = geo_data.global_qbx_centers() icontaining_tgt_box_vec = qbx_center_to_target_box[tgt_icenter_vec] @@ -391,8 +389,6 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): del tgt_icenter del icontaining_tgt_box - print("end par data prep") - if self.dim == 3: # This gets max'd onto: pass initialized version. ier = np.zeros(ngqbx_centers, dtype=np.int32) -- GitLab From beb261f94094ae1d8d5bb0c00cac55ff7ff3332d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 23:02:26 -0500 Subject: [PATCH 22/53] Remove curl_S --- pytential/symbolic/pde/maxwell/__init__.py | 18 +++++++++++++----- pytential/symbolic/primitives.py | 4 ---- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 93dfffc6..bc47a53e 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -49,7 +49,9 @@ class PECAugmentedMFIEOperator: def j_operator(self, loc, Jt): Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") return xyz_to_tangential( - (loc*0.5)*Jxyz - sym.nxcurl_S(self.kernel, 0, Jxyz, k=self.k)) + (loc*0.5)*Jxyz - sym.n_cross( + sym.curl(sym.S(self.kernel, Jxyz, k=self.k, + qbx_forced_limit="avg")))) def j_rhs(self, Hinc_xyz): return xyz_to_tangential(sym.n_cross(Hinc_xyz)) @@ -72,7 +74,7 @@ class PECAugmentedMFIEOperator: # use - n x n x v = v_tangential 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 + H_scat = sym.curl(sym.S(self.kernel, Jxyz, k=self.k)) + (loc*0.5)*Jxyz return sym.join_fields(E_scat, H_scat) @@ -83,7 +85,7 @@ class PECAugmentedMFIEOperator: grad_phi = sym.grad(3, sym.S(self.kernel, rho, k=self.k)) E_scat = 1j*self.k*A - grad_phi - H_scat = sym.curl_S(self.kernel, Jxyz, k=self.k) + H_scat = sym.curl(sym.S(self.kernel, Jxyz, k=self.k)) return sym.join_fields(E_scat, H_scat) @@ -131,7 +133,10 @@ class MuellerAugmentedMFIEOperator(object): k0, k1 = self.ks S = partial(sym.S, self.kernel, qbx_forced_limit="avg") - curl_S = partial(sym.curl_S, self.kernel, qbx_forced_limit="avg") + + def curl_S(dens): + return sym.curl(sym.S(self.kernel, dens, qbx_forced_limit="avg")) + grad = partial(sym.grad, 3) E0 = sym.cse(1j*omega*mu0*eps0*S(Jxyz, k=k0) + @@ -173,7 +178,10 @@ 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=None, k=k) + + def curl_S(dens): + return sym.curl(sym.S(self.kernel, dens, qbx_forced_limit=None, k=k)) + grad = partial(sym.grad, 3) E0 = 1j*k*eps*S(Jxyz) + mu*curl_S(Mxyz) - grad(S(u.rho_e)) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 12c96d8c..4862d57e 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -1044,10 +1044,6 @@ def curl(vec): for m in range(3) for n in range(3)) for l in range(3)]) - -def curl_S(kernel, arg, **kwargs): # noqa: N802 - return curl(S(kernel, arg, **kwargs)) - # }}} -- GitLab From 171312fb2bfeaceb38605ae117926a7e32fc1a93 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 23:02:54 -0500 Subject: [PATCH 23/53] Maxwell tests executes, does not solve PDE --- test/test_maxwell.py | 155 ++++++++++++++++++++++++++++--------------- 1 file changed, 103 insertions(+), 52 deletions(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index f652ae65..0e74c570 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -33,6 +33,7 @@ from pytential import bind, sym from sumpy.visualization import FieldPlotter # noqa from sumpy.point_calculus import CalculusPatch, frequency_domain_maxwell +from sumpy.tools import vector_from_device from pytential.target import PointsTarget import logging @@ -43,6 +44,8 @@ from pyopencl.tools import ( # noqa class TestCase: + fmm_backend = "fmmlib" + def __init__(self, k, is_interior, resolutions, qbx_order, fmm_order): self.k = k @@ -53,6 +56,9 @@ class TestCase: class SphereTestCase(TestCase): + target_order = 8 + gmres_tol = 1e-10 + def get_mesh(self, resolution, target_order): from meshmode.mesh.generation import generate_icosphere from meshmode.mesh.refinement import refine_uniformly @@ -82,9 +88,9 @@ class SphereTestCase(TestCase): cl.array.to_device(queue, sources)) -tc_int = SphereTestCase(k=10.2, is_interior=True, resolutions=[0, 1], +tc_int = SphereTestCase(k=1.2, is_interior=True, resolutions=[0, 1], qbx_order=3, fmm_order=10) -tc_ext = SphereTestCase(k=10.2, is_interior=False, resolutions=[0, 1], +tc_ext = SphereTestCase(k=1.2, is_interior=False, resolutions=[0, 1], qbx_order=3, fmm_order=10) @@ -105,58 +111,55 @@ def get_sym_maxwell_source(kernel, jxyz, k): tc_ext ]) def test_mfie_from_source(ctx_getter, case, visualize=False): - logging.basicConfig(level=logging.INFO) + #logging.basicConfig(level=logging.INFO) cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) pytest.importorskip("pyfmmlib") - pytest.importorskip("scipy") np.random.seed(12) + knl_kwargs = {"k": case.k} + # {{{ come up with a solution to Maxwell's equations j_sym = sym.make_sym_vector("j", 3) + jt_sym = sym.make_sym_vector("jt", 2) + rho_sym = sym.var("rho") from pytential.symbolic.pde.maxwell import PECAugmentedMFIEOperator mfie = PECAugmentedMFIEOperator() test_source = case.get_source(queue) - calc_patch = CalculusPatch(np.array([-2, 0, 0]), h=0.01) + calc_patch = CalculusPatch(np.array([-3, 0, 0]), h=0.01) calc_patch_tgt = PointsTarget(calc_patch.points) rng = cl.clrandom.PhiloxGenerator(cl_ctx, seed=12) src_j = rng.normal(queue, (3, test_source.nnodes), dtype=np.float64) - pde_test_field = bind( - (test_source, calc_patch_tgt), - get_sym_maxwell_source(mfie.kernel, j_sym, mfie.k) - )(queue, j=src_j, k=case.k) + def eval_inc_field_at(tgt): + return bind( + (test_source, tgt), + get_sym_maxwell_source(mfie.kernel, j_sym, mfie.k) + )(queue, j=src_j, k=case.k) - from sumpy.tools import vector_from_device - pde_test_e = vector_from_device(queue, pde_test_field[:3]) - pde_test_h = vector_from_device(queue, pde_test_field[3:]) + pde_test_inc = eval_inc_field_at(calc_patch_tgt) - print([ + pde_test_e = vector_from_device(queue, pde_test_inc[:3]) + pde_test_h = vector_from_device(queue, pde_test_inc[3:]) + + assert max( calc_patch.norm(x, np.inf) / calc_patch.norm(pde_test_e, np.inf) for x in frequency_domain_maxwell( - calc_patch, pde_test_e, pde_test_h, case.k)]) - - - 1/0 + calc_patch, pde_test_e, pde_test_h, case.k)) < 1e-6 + loc_sign = -1 if case.is_interior else +1 + # xyz_mfie_bdry_vol_op = mfie.boundary_field(0, Jxyz) + # xyz_mfie_vol_op = mfie.volume_field(Jxyz) - - - - mfie_ext_op = xyz_to_tangential(mfie.boundary_field(+1, Jxyz)) - mfie_int_op = xyz_to_tangential(mfie.boundary_field(-1, Jxyz)) - xyz_mfie_bdry_vol_op = mfie.boundary_field(0, Jxyz) - xyz_mfie_vol_op = mfie.volume_field(Jxyz) - - n_cross_op = n_cross(make_vector_field("vec", 3)) + # n_cross_op = n_cross(make_vector_field("vec", 3)) from pytools.convergence import EOCRecorder @@ -164,46 +167,94 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): eoc_rec_E = EOCRecorder() eoc_rec_H = EOCRecorder() + from pytential.qbx import QBXLayerPotentialSource + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory for resolution in case.resolutions: - scat_mesh = generate_icosphere(2, i) - if is_interior: - tgt_mesh = generate_icosphere(5, 0) - else: - tgt_mesh = generate_icosphere(0.5, 1) + scat_mesh = case.get_mesh(resolution, case.target_order) + test_mesh = case.get_test_mesh(case.target_order) + + pre_scat_discr = Discretization( + cl_ctx, scat_mesh, + InterpolatoryQuadratureSimplexGroupFactory(case.target_order)) + qbx, _ = QBXLayerPotentialSource( + pre_scat_discr, fine_order=4*case.target_order, + qbx_order=case.qbx_order, + fmm_order=case.fmm_order, fmm_backend=case.fmm_backend + ).with_refinement() + scat_discr = qbx.density_discr + test_discr = Discretization( + cl_ctx, test_mesh, + InterpolatoryQuadratureSimplexGroupFactory(case.target_order)) + + inc_field_scat = eval_inc_field_at(scat_discr) + inc_field_test = eval_inc_field_at(test_discr) + + # nxHinc_xyz_mid = sbind(n_cross_op)(vec=H_scat) + # nxH_tgt_true = bind(n_cross_op, tgt_discr, iprec=iprec)( + # vec=H_tgt_true) - scat_discr = FlatConstantDiscretization(scat_mesh) - tgt_discr = FlatConstantDiscretization(tgt_mesh) + # {{{ system solve - def sbind(op): return bind(op, scat_discr, iprec=iprec) + # nxHinc_t_mid = sbind( + # xyz_to_tangential(make_vector_field("vec", 3)) + # )(vec=nxHinc_xyz_mid) - source_func = gen_field_multipole + inc_xyz_sym = sym.make_sym_vector("inc_fld", 6) + e_inc_xyz_sym = inc_xyz_sym[:3] + h_inc_xyz_sym = inc_xyz_sym[3:] - E_scat, H_scat = source_func(k, is_interior, scat_mesh.centroids) - E_tgt_true, H_tgt_true = source_func(k, is_interior, tgt_mesh.centroids) + bound_j_op = bind(qbx, mfie.j_operator(loc_sign, jt_sym)) + j_rhs = bind(qbx, mfie.j_rhs(h_inc_xyz_sym))( + queue, inc_fld=inc_field_scat, **knl_kwargs) - nxHinc_xyz_mid = sbind(n_cross_op)(vec=H_scat) - nxH_tgt_true = bind(n_cross_op, tgt_discr, iprec=iprec)( - vec=H_tgt_true) + from pytential.solve import gmres + gmres_result = gmres( + bound_j_op.scipy_op(queue, "jt", np.complex128, **knl_kwargs), + j_rhs, + tol=case.gmres_tol, + progress=True, + hard_failure=True) - # {{{ system solve + jt = gmres_result.solution - nxHinc_t_mid = sbind( - xyz_to_tangential(make_vector_field("vec", 3)) - )(vec=nxHinc_xyz_mid) + bound_rho_op = bind(qbx, mfie.rho_operator(loc_sign, rho_sym)) + rho_rhs = bind(qbx, mfie.rho_rhs(jt_sym, e_inc_xyz_sym))( + queue, jt=jt, inc_fld=inc_field_scat, **knl_kwargs) - if is_interior: - mfie_solve_op = mfie_int_op - else: - mfie_solve_op = mfie_ext_op + gmres_result = gmres( + bound_rho_op.scipy_op(queue, "rho", np.complex128, **knl_kwargs), + rho_rhs, + tol=case.gmres_tol, + progress=True, + hard_failure=True) - Jt_mid = solve_lin_op( - sbind(mfie_solve_op).scipy_op("Jt"), - nxHinc_t_mid) + rho = gmres_result.solution # }}} - # {{{ error in nxH + # {{{ volume eval + + sym_repr = mfie.scattered_volume_field(jt_sym, rho_sym) + + def eval_repr_at(tgt): + return bind((qbx, tgt), sym_repr)(queue, jt=jt, rho=rho, **knl_kwargs) + + pde_test_repr = eval_repr_at(calc_patch_tgt) + + pde_test_e = vector_from_device(queue, pde_test_repr[:3]) + pde_test_h = vector_from_device(queue, pde_test_repr[3:]) + + print([ + calc_patch.norm(x, np.inf) / calc_patch.norm(pde_test_e, np.inf) + for x in frequency_domain_maxwell( + calc_patch, pde_test_e, pde_test_h, case.k)]) + + calc_patch_tgt + + nxH_tgt = bind(xyz_mfie_bdry_vol_op, (scat_discr, tgt_discr), iprec=iprec) \ (Jt=Jt_mid) -- GitLab From c031541a4748b567755f3b83e68e9784aa11b11c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 23:11:05 -0500 Subject: [PATCH 24/53] MFIE representation sign fixes --- pytential/symbolic/pde/maxwell/__init__.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index bc47a53e..7e9263d4 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -63,7 +63,7 @@ class PECAugmentedMFIEOperator: Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") return (sym.n_dot(Einc_xyz) - + 1j*self.k*sym.n_dot(sym.S(self.kernel, Jxyz, k=self.k))) + - 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") @@ -73,7 +73,7 @@ class PECAugmentedMFIEOperator: # use - n x n x v = v_tangential - E_scat = 1j*self.k*A - grad_phi + 0.5*loc*rho + E_scat = - 1j*self.k*A - grad_phi + 0.5*loc*rho H_scat = sym.curl(sym.S(self.kernel, Jxyz, k=self.k)) + (loc*0.5)*Jxyz return sym.join_fields(E_scat, H_scat) @@ -84,7 +84,7 @@ class PECAugmentedMFIEOperator: A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=None) grad_phi = sym.grad(3, sym.S(self.kernel, rho, k=self.k)) - E_scat = 1j*self.k*A - grad_phi + E_scat = - 1j*self.k*A - grad_phi H_scat = sym.curl(sym.S(self.kernel, Jxyz, k=self.k)) return sym.join_fields(E_scat, H_scat) -- GitLab From 0418c5e72540086a00071bff8681a841f889052f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 5 Sep 2017 23:59:38 -0500 Subject: [PATCH 25/53] More work towards MFIE test --- test/test_maxwell.py | 65 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 62 insertions(+), 3 deletions(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index 0e74c570..d63b3fcd 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -235,12 +235,17 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): # }}} + jxyz = bind(qbx, sym.tangential_to_xyz(jt_sym))(queue, jt=jt) + # {{{ volume eval sym_repr = mfie.scattered_volume_field(jt_sym, rho_sym) - def eval_repr_at(tgt): - return bind((qbx, tgt), sym_repr)(queue, jt=jt, rho=rho, **knl_kwargs) + def eval_repr_at(tgt, source=None): + if source is None: + source = qbx + + return bind((source, tgt), sym_repr)(queue, jt=jt, rho=rho, **knl_kwargs) pde_test_repr = eval_repr_at(calc_patch_tgt) @@ -252,7 +257,61 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): for x in frequency_domain_maxwell( calc_patch, pde_test_e, pde_test_h, case.k)]) - calc_patch_tgt + #test_repr = eval_repr_at(test_discr) + + # {{{ visualization + + if visualize: + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, scat_discr, case.target_order+3) + + bdry_normals = bind(scat_discr, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ + ("j", jxyz), + ("rho", rho), + ("Einc", inc_field_scat[:3]), + ("Hinc", inc_field_scat[3:]), + ("bdry_normals", bdry_normals), + ]) + + bbox_center = np.zeros(3) + bbox_size = 6 + + fplot = FieldPlotter( + bbox_center, extent=bbox_size, npoints=(150, 150, 5)) + + from pytential.qbx import QBXTargetAssociationFailedException + + qbx_tgt_tol = qbx.copy(target_association_tolerance=0.15) + + fplot_tgt = PointsTarget(fplot.points) + try: + fplot_repr = eval_repr_at(fplot_tgt, source=qbx_tgt_tol) + except QBXTargetAssociationFailedException as e: + fplot.write_vtk_file( + "failed-targets.vts", + [ + ("failed_targets", e.failed_target_flags.get(queue)) + ]) + raise + + fplot_inc = eval_inc_field_at(fplot_tgt) + + fplot.write_vtk_file( + "potential.vts", + [ + ("E", fplot_repr[:3]), + ("H", fplot_repr[3:]), + ("Einc", fplot_inc[:3]), + ("Hinc", fplot_inc[3:]), + ] + ) + + # }}} + + -- GitLab From ffc32497d06a6724f5fce57571b427300d2651a1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 6 Sep 2017 00:01:53 -0500 Subject: [PATCH 26/53] Fix MFIE notes --- contrib/notes/mfie.tm | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/contrib/notes/mfie.tm b/contrib/notes/mfie.tm index ff3d86b7..b2f2e648 100644 --- a/contrib/notes/mfie.tm +++ b/contrib/notes/mfie.tm @@ -1,4 +1,4 @@ - + @@ -65,26 +65,29 @@ n\+E|)>=\. - Now use the representation =i*k*A-\\=i*k*SJ-\S\> + Now use the representation =-i*k*A-\\=i*k*SJ-\S\> and obtain <\equation*> - n\E+J-\S\|)>|)>=\. + n\E+J-\S\|)>|)>=\. Carrying out the limit, we obtain: <\equation*> - n\E+n\J|)>-S\+\=\. + n\E-n\J|)>-S\+\=\. Rearrange: <\equation*> - n\E+n\J|)>=\+S\. + n\E-n\J|)>=\+S\. + +> + <\references> <\collection> > @@ -96,15 +99,15 @@ <\auxiliary> <\collection> <\associate|toc> - |math-font-series||Augmented + |math-font-series||1Augmented MFIE derivation> |.>>>>|> - |MFIE derivation + |1.1MFIE derivation |.>>>>|> > - ||\> + |1.2|\> postprocessor derivation |.>>>>|> > -- GitLab From b9a64594d5a338e11da5c6d1de5d6ed9303f24a5 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 6 Sep 2017 00:31:05 -0500 Subject: [PATCH 27/53] Promote unspecified qbx_forced_limit to a proper warning, from a DeprecationWarning --- pytential/symbolic/primitives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 4862d57e..d00aef68 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -847,7 +847,7 @@ def S(kernel, density, # noqa if qbx_forced_limit is _unspecified: warn("not specifying qbx_forced_limit on call to 'S' is deprecated, " - "defaulting to +1", DeprecationWarning, stacklevel=2) + "defaulting to +1", stacklevel=2) qbx_forced_limit = +1 return IntG(kernel, density, qbx_forced_limit, source, target, -- GitLab From 90de8ffcb7e863c900f6247b2a19bee243a2600a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 6 Sep 2017 00:31:37 -0500 Subject: [PATCH 28/53] MFIE representation: Use QBX vol eval --- pytential/symbolic/pde/maxwell/__init__.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 7e9263d4..1ef6b991 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -82,10 +82,12 @@ class PECAugmentedMFIEOperator: Jxyz = sym.cse(sym.tangential_to_xyz(Jt), "Jxyz") A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=None) - grad_phi = sym.grad(3, sym.S(self.kernel, rho, k=self.k)) + grad_phi = sym.grad(3, sym.S(self.kernel, rho, k=self.k, + qbx_forced_limit=None)) E_scat = - 1j*self.k*A - grad_phi - H_scat = sym.curl(sym.S(self.kernel, Jxyz, k=self.k)) + H_scat = sym.curl(sym.S(self.kernel, Jxyz, k=self.k, + qbx_forced_limit=None)) return sym.join_fields(E_scat, H_scat) -- GitLab From f0b2d7cc06a206e65b75904d30646ea6fc8e1677 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 15 Sep 2017 00:01:09 -0500 Subject: [PATCH 29/53] Initial version of MFIE showing E/H convergence (but: BC checks etc remain) --- test/test_maxwell.py | 111 ++++++++++++++++--------------------------- 1 file changed, 42 insertions(+), 69 deletions(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index d63b3fcd..a530bfea 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -29,7 +29,7 @@ import pyopencl.clmath # noqa import pyopencl.clrandom # noqa import pytest -from pytential import bind, sym +from pytential import bind, sym, norm from sumpy.visualization import FieldPlotter # noqa from sumpy.point_calculus import CalculusPatch, frequency_domain_maxwell @@ -66,7 +66,7 @@ class SphereTestCase(TestCase): generate_icosphere(2, target_order), resolution) - def get_test_mesh(self, target_order): + def get_observation_mesh(self, target_order): from meshmode.mesh.generation import generate_icosphere if self.is_interior: @@ -81,7 +81,9 @@ class SphereTestCase(TestCase): else: source_ctr = np.array([[5, 0.1, 0.15]]).T - sources = source_ctr + np.random.rand(3, 10) + source_rad = 0.3 + + sources = source_ctr + source_rad*2*(np.random.rand(3, 10)-0.5) from pytential.source import PointPotentialSource return PointPotentialSource( queue.context, @@ -111,7 +113,7 @@ def get_sym_maxwell_source(kernel, jxyz, k): tc_ext ]) def test_mfie_from_source(ctx_getter, case, visualize=False): - #logging.basicConfig(level=logging.INFO) + logging.basicConfig(level=logging.INFO) cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) @@ -164,8 +166,8 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): from pytools.convergence import EOCRecorder eoc_rec_nxH = EOCRecorder() - eoc_rec_E = EOCRecorder() - eoc_rec_H = EOCRecorder() + eoc_rec_e = EOCRecorder() + eoc_rec_h = EOCRecorder() from pytential.qbx import QBXLayerPotentialSource from meshmode.discretization import Discretization @@ -174,7 +176,7 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): for resolution in case.resolutions: scat_mesh = case.get_mesh(resolution, case.target_order) - test_mesh = case.get_test_mesh(case.target_order) + observation_mesh = case.get_observation_mesh(case.target_order) pre_scat_discr = Discretization( cl_ctx, scat_mesh, @@ -185,12 +187,12 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): fmm_order=case.fmm_order, fmm_backend=case.fmm_backend ).with_refinement() scat_discr = qbx.density_discr - test_discr = Discretization( - cl_ctx, test_mesh, + obs_discr = Discretization( + cl_ctx, observation_mesh, InterpolatoryQuadratureSimplexGroupFactory(case.target_order)) inc_field_scat = eval_inc_field_at(scat_discr) - inc_field_test = eval_inc_field_at(test_discr) + inc_field_obs = eval_inc_field_at(obs_discr) # nxHinc_xyz_mid = sbind(n_cross_op)(vec=H_scat) # nxH_tgt_true = bind(n_cross_op, tgt_discr, iprec=iprec)( @@ -257,7 +259,7 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): for x in frequency_domain_maxwell( calc_patch, pde_test_e, pde_test_h, case.k)]) - #test_repr = eval_repr_at(test_discr) + # }}} # {{{ visualization @@ -284,7 +286,7 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): from pytential.qbx import QBXTargetAssociationFailedException - qbx_tgt_tol = qbx.copy(target_association_tolerance=0.15) + qbx_tgt_tol = qbx.copy(target_association_tolerance=0.2) fplot_tgt = PointsTarget(fplot.points) try: @@ -302,87 +304,58 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): fplot.write_vtk_file( "potential.vts", [ - ("E", fplot_repr[:3]), - ("H", fplot_repr[3:]), - ("Einc", fplot_inc[:3]), - ("Hinc", fplot_inc[3:]), + ("E", vector_from_device(queue, fplot_repr[:3])), + ("H", vector_from_device(queue, fplot_repr[3:])), + ("Einc", vector_from_device(queue, fplot_inc[:3])), + ("Hinc", vector_from_device(queue, fplot_inc[3:])), ] ) # }}} + # {{{ error in E, H + obs_repr = eval_repr_at(obs_discr) + obs_e = obs_repr[:3] + obs_h = obs_repr[3:] + inc_obs_e = inc_field_obs[:3] + inc_obs_h = inc_field_obs[3:] - nxH_tgt = bind(xyz_mfie_bdry_vol_op, (scat_discr, tgt_discr), iprec=iprec) \ - (Jt=Jt_mid) - rel_err_nxH = (tgt_discr.norm(nxH_tgt_true+nxH_tgt) - / tgt_discr.norm(nxH_tgt_true)) - - # }}} - - # {{{ error in E, H + def obs_norm(f): + return norm(obs_discr, queue, f, p=np.inf) - EH_tgt = bind(xyz_mfie_vol_op, (scat_discr, tgt_discr), iprec=iprec) \ - (Jt=Jt_mid) - E_tgt = EH_tgt[:3] - H_tgt = EH_tgt[3:] - rel_err_H = (tgt_discr.norm(H_tgt_true-H_tgt) - / tgt_discr.norm(H_tgt_true)) - rel_err_E = (tgt_discr.norm(E_tgt_true-E_tgt) - / tgt_discr.norm(E_tgt_true)) + # FIXME: Why "+"? + rel_err_e = (obs_norm(inc_obs_e + obs_e) + / obs_norm(inc_obs_e)) + rel_err_h = (obs_norm(inc_obs_h + obs_h) + / obs_norm(inc_obs_h)) # }}} - print("ERR", h, rel_err_nxH, rel_err_H, rel_err_E) + h_max = qbx.h_max + print("ERR", h_max, rel_err_h, rel_err_e) - eoc_rec_nxH.add_data_point(h, rel_err_nxH) - eoc_rec_H.add_data_point(h, rel_err_H) - eoc_rec_E.add_data_point(h, rel_err_E) + eoc_rec_h.add_data_point(h_max, rel_err_h) + eoc_rec_e.add_data_point(h_max, rel_err_e) - # {{{ visualization - if write_vis: - from pyvisfile.silo import SiloFile - from hellskitchen.visualization import add_to_silo_file - from pytools.obj_array import oarray_real_copy - - Jxyz_mid = bind( - tangential_to_xyz(make_vector_field("vec", 2)))(vec=Jt_mid) - - silo = SiloFile("mfie-int%s-%d.silo" % (is_interior, i)) - add_to_silo_file(silo, scat_geo, cell_data=[ - ("nxHinc_mid", oarray_real_copy(nxHinc_xyz_mid)), - ("J_mid", oarray_real_copy(Jxyz_mid)), - ]) - add_to_silo_file(silo, tgt_geo, cell_data=[ - ("nxHinc_tgt", oarray_real_copy(nxH_tgt)), - ("nxHinc_tgt_true", oarray_real_copy(nxH_tgt_true)), - ("Hinc_tgt", oarray_real_copy(H_tgt)), - ("Hinc_tgt_true", oarray_real_copy(H_tgt_true)), - ("Einc_tgt", oarray_real_copy(E_tgt)), - ("Einc_tgt_true", oarray_real_copy(E_tgt_true)), - ], mesh_name="tgt_mesh") - silo.close() - - # }}} + # TODO: Add Maxwellian-ness convergence + # TODO: Check that total field verifies BC + # TODO: Check for extinction on the interior print("--------------------------------------------------------") print("is_interior=%s" % case.is_interior) print("--------------------------------------------------------") for which_eoc, eoc_rec in [ - ("nxH", eoc_rec_nxH), - ("H", eoc_rec_H), - ("E", eoc_rec_E)]: + #("nxH", eoc_rec_nxH), + ("H", eoc_rec_h), + ("E", eoc_rec_e)]: print(which_eoc) print(eoc_rec.pretty_print()) if len(eoc_rec.history) > 1: - assert eoc_rec.order_estimate() > 0.8 - - - - + assert eoc_rec.order_estimate() > case.qbx_order - 1 # You can test individual routines by typing -- GitLab From 960d4e82874b3098e58756163b448dd78a1663e1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 16 Sep 2017 09:31:44 -0500 Subject: [PATCH 30/53] MFIE test: Add EOC for Maxwellian-ness --- test/test_maxwell.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index a530bfea..b5afc6f5 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -165,7 +165,7 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): from pytools.convergence import EOCRecorder - eoc_rec_nxH = EOCRecorder() + eoc_rec_repr_maxwell = EOCRecorder() eoc_rec_e = EOCRecorder() eoc_rec_h = EOCRecorder() @@ -186,6 +186,8 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): qbx_order=case.qbx_order, fmm_order=case.fmm_order, fmm_backend=case.fmm_backend ).with_refinement() + h_max = qbx.h_max + scat_discr = qbx.density_discr obs_discr = Discretization( cl_ctx, observation_mesh, @@ -254,10 +256,13 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): pde_test_e = vector_from_device(queue, pde_test_repr[:3]) pde_test_h = vector_from_device(queue, pde_test_repr[3:]) - print([ + maxwell_residuals = [ calc_patch.norm(x, np.inf) / calc_patch.norm(pde_test_e, np.inf) for x in frequency_domain_maxwell( - calc_patch, pde_test_e, pde_test_h, case.k)]) + calc_patch, pde_test_e, pde_test_h, case.k)] + print("Maxwell residuals:", maxwell_residuals) + + eoc_rec_repr_maxwell.add_data_point(h_max, max(maxwell_residuals)) # }}} @@ -334,28 +339,31 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): # }}} - h_max = qbx.h_max print("ERR", h_max, rel_err_h, rel_err_e) eoc_rec_h.add_data_point(h_max, rel_err_h) eoc_rec_e.add_data_point(h_max, rel_err_e) - # TODO: Add Maxwellian-ness convergence # TODO: Check that total field verifies BC # TODO: Check for extinction on the interior print("--------------------------------------------------------") print("is_interior=%s" % case.is_interior) print("--------------------------------------------------------") + + good = True for which_eoc, eoc_rec in [ - #("nxH", eoc_rec_nxH), + ("maxwell", eoc_rec_repr_maxwell), ("H", eoc_rec_h), ("E", eoc_rec_e)]: print(which_eoc) print(eoc_rec.pretty_print()) if len(eoc_rec.history) > 1: - assert eoc_rec.order_estimate() > case.qbx_order - 1 + if eoc_rec.order_estimate() < case.qbx_order - 1: + good = False + + assert good # You can test individual routines by typing -- GitLab From a14854b5e6d5234d5ec85b9734b5a69bbd545c24 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 16 Sep 2017 09:42:54 -0500 Subject: [PATCH 31/53] MFIE test: Adjust EOC tol for Maxwellian-ness --- test/test_maxwell.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index b5afc6f5..5a5bb14d 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -157,6 +157,8 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): for x in frequency_domain_maxwell( calc_patch, pde_test_e, pde_test_h, case.k)) < 1e-6 + # }}} + loc_sign = -1 if case.is_interior else +1 # xyz_mfie_bdry_vol_op = mfie.boundary_field(0, Jxyz) # xyz_mfie_vol_op = mfie.volume_field(Jxyz) @@ -352,15 +354,15 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): print("--------------------------------------------------------") good = True - for which_eoc, eoc_rec in [ - ("maxwell", eoc_rec_repr_maxwell), - ("H", eoc_rec_h), - ("E", eoc_rec_e)]: + for which_eoc, eoc_rec, order_tol in [ + ("maxwell", eoc_rec_repr_maxwell, 1.5), + ("H", eoc_rec_h, 1), + ("E", eoc_rec_e, 1)]: print(which_eoc) print(eoc_rec.pretty_print()) if len(eoc_rec.history) > 1: - if eoc_rec.order_estimate() < case.qbx_order - 1: + if eoc_rec.order_estimate() < case.qbx_order - order_tol: good = False assert good -- GitLab From 202e24e3406cd4080b6588ed2a10cbee6e4d4ff9 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 16 Sep 2017 09:44:34 -0500 Subject: [PATCH 32/53] MFIE test vis: use unique file names --- test/test_maxwell.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index 5a5bb14d..8d72d8a4 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -309,7 +309,7 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): fplot_inc = eval_inc_field_at(fplot_tgt) fplot.write_vtk_file( - "potential.vts", + "potential-%s.vts" % resolution, [ ("E", vector_from_device(queue, fplot_repr[:3])), ("H", vector_from_device(queue, fplot_repr[3:])), -- GitLab From 47637a8b04beb34d24b14a4e145037cf672942dd Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 16 Sep 2017 11:48:42 -0500 Subject: [PATCH 33/53] MFIE: Adjust to exp(-i omega t) time dep convention --- pytential/symbolic/pde/maxwell/__init__.py | 14 ++++++-------- test/test_maxwell.py | 10 +++++++--- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 1ef6b991..0214aeb1 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -63,17 +63,17 @@ class PECAugmentedMFIEOperator: Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") return (sym.n_dot(Einc_xyz) - - 1j*self.k*sym.n_dot(sym.S(self.kernel, Jxyz, k=self.k))) + + 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 = sym.S(self.kernel, Jxyz, k=self.k) - grad_phi = sym.grad(3, sym.S(self.kernel, rho, k=self.k)) + phi = sym.S(self.kernel, rho, k=self.k) # use - n x n x v = v_tangential - E_scat = - 1j*self.k*A - grad_phi + 0.5*loc*rho + E_scat = 1j*self.k*A - sym.grad(3, phi) + 0.5*loc*rho H_scat = sym.curl(sym.S(self.kernel, Jxyz, k=self.k)) + (loc*0.5)*Jxyz return sym.join_fields(E_scat, H_scat) @@ -82,12 +82,10 @@ class PECAugmentedMFIEOperator: Jxyz = sym.cse(sym.tangential_to_xyz(Jt), "Jxyz") A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=None) - grad_phi = sym.grad(3, sym.S(self.kernel, rho, k=self.k, - qbx_forced_limit=None)) + phi = sym.S(self.kernel, rho, k=self.k, qbx_forced_limit=None) - E_scat = - 1j*self.k*A - grad_phi - H_scat = sym.curl(sym.S(self.kernel, Jxyz, k=self.k, - qbx_forced_limit=None)) + E_scat = 1j*self.k*A - sym.grad(3, phi) + H_scat = sym.curl(A) return sym.join_fields(E_scat, H_scat) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index 8d72d8a4..635c93b6 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -103,8 +103,10 @@ def get_sym_maxwell_source(kernel, jxyz, k): A = sym.curl(sym.S(kernel, jxyz, k=k, qbx_forced_limit=None)) # https://en.wikipedia.org/w/index.php?title=Maxwell%27s_equations&oldid=798940325#Alternative_formulations + # (Vector calculus/Potentials/Any Gauge) + # assumed time dependence exp(-1j*omega*t) return sym.join_fields( - - 1j*k*A, + 1j*k*A, sym.curl(A)) @@ -152,10 +154,12 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): pde_test_e = vector_from_device(queue, pde_test_inc[:3]) pde_test_h = vector_from_device(queue, pde_test_inc[3:]) - assert max( + source_maxwell_resids = [ calc_patch.norm(x, np.inf) / calc_patch.norm(pde_test_e, np.inf) for x in frequency_domain_maxwell( - calc_patch, pde_test_e, pde_test_h, case.k)) < 1e-6 + calc_patch, pde_test_e, pde_test_h, case.k)] + print("Source Maxwell residuals:", source_maxwell_resids) + assert max(source_maxwell_resids) < 1e-6 # }}} -- GitLab From 2d8fe0ed5cffad769ae96cdea62b0bbab67c3bf1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 16 Sep 2017 11:53:26 -0500 Subject: [PATCH 34/53] MFIE op: Add comment about shape of time dependency. --- pytential/symbolic/pde/maxwell/__init__.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 0214aeb1..030c4185 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -35,10 +35,13 @@ cse = sym.cse # {{{ Charge-Current MFIE class PECAugmentedMFIEOperator: - """Magnetic Field Integral Equation operator, + r"""Magnetic Field Integral Equation operator, under the assumption of no surface charges. see :file:`contrib/notes/mfie.tm` + + Uses the sign convention :math:`\exp(-1 \omega t)` + for the sinusoidal time dependency. """ def __init__(self, k=sym.var("k")): -- GitLab From 10c2bd4312f00756f9c2d0c76f90556279eaccbe Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 16 Sep 2017 17:30:44 -0500 Subject: [PATCH 35/53] Use EHField object to shorten MFIE code some --- test/test_maxwell.py | 92 +++++++++++++++++++++----------------------- 1 file changed, 44 insertions(+), 48 deletions(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index 635c93b6..26ec2e70 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -110,6 +110,20 @@ def get_sym_maxwell_source(kernel, jxyz, k): sym.curl(A)) +class EHField(object): + def __init__(self, eh_field): + assert len(eh_field) == 6 + self.field = eh_field + + @property + def e(self): + return self.field[:3] + + @property + def h(self): + return self.field[3:] + + @pytest.mark.parametrize("case", [ tc_int, tc_ext @@ -149,15 +163,13 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): get_sym_maxwell_source(mfie.kernel, j_sym, mfie.k) )(queue, j=src_j, k=case.k) - pde_test_inc = eval_inc_field_at(calc_patch_tgt) - - pde_test_e = vector_from_device(queue, pde_test_inc[:3]) - pde_test_h = vector_from_device(queue, pde_test_inc[3:]) + pde_test_inc = EHField( + vector_from_device(queue, eval_inc_field_at(calc_patch_tgt))) source_maxwell_resids = [ - calc_patch.norm(x, np.inf) / calc_patch.norm(pde_test_e, np.inf) + calc_patch.norm(x, np.inf) / calc_patch.norm(pde_test_inc.e, np.inf) for x in frequency_domain_maxwell( - calc_patch, pde_test_e, pde_test_h, case.k)] + calc_patch, pde_test_inc.e, pde_test_inc.h, case.k)] print("Source Maxwell residuals:", source_maxwell_resids) assert max(source_maxwell_resids) < 1e-6 @@ -199,26 +211,16 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): cl_ctx, observation_mesh, InterpolatoryQuadratureSimplexGroupFactory(case.target_order)) - inc_field_scat = eval_inc_field_at(scat_discr) - inc_field_obs = eval_inc_field_at(obs_discr) - - # nxHinc_xyz_mid = sbind(n_cross_op)(vec=H_scat) - # nxH_tgt_true = bind(n_cross_op, tgt_discr, iprec=iprec)( - # vec=H_tgt_true) + inc_field_scat = EHField(eval_inc_field_at(scat_discr)) + inc_field_obs = EHField(eval_inc_field_at(obs_discr)) # {{{ system solve - # nxHinc_t_mid = sbind( - # xyz_to_tangential(make_vector_field("vec", 3)) - # )(vec=nxHinc_xyz_mid) - - inc_xyz_sym = sym.make_sym_vector("inc_fld", 6) - e_inc_xyz_sym = inc_xyz_sym[:3] - h_inc_xyz_sym = inc_xyz_sym[3:] + inc_xyz_sym = EHField(sym.make_sym_vector("inc_fld", 6)) bound_j_op = bind(qbx, mfie.j_operator(loc_sign, jt_sym)) - j_rhs = bind(qbx, mfie.j_rhs(h_inc_xyz_sym))( - queue, inc_fld=inc_field_scat, **knl_kwargs) + j_rhs = bind(qbx, mfie.j_rhs(inc_xyz_sym.h))( + queue, inc_fld=inc_field_scat.field, **knl_kwargs) from pytential.solve import gmres gmres_result = gmres( @@ -231,8 +233,8 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): jt = gmres_result.solution bound_rho_op = bind(qbx, mfie.rho_operator(loc_sign, rho_sym)) - rho_rhs = bind(qbx, mfie.rho_rhs(jt_sym, e_inc_xyz_sym))( - queue, jt=jt, inc_fld=inc_field_scat, **knl_kwargs) + rho_rhs = bind(qbx, mfie.rho_rhs(jt_sym, inc_xyz_sym.e))( + queue, jt=jt, inc_fld=inc_field_scat.field, **knl_kwargs) gmres_result = gmres( bound_rho_op.scipy_op(queue, "rho", np.complex128, **knl_kwargs), @@ -257,15 +259,13 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): return bind((source, tgt), sym_repr)(queue, jt=jt, rho=rho, **knl_kwargs) - pde_test_repr = eval_repr_at(calc_patch_tgt) - - pde_test_e = vector_from_device(queue, pde_test_repr[:3]) - pde_test_h = vector_from_device(queue, pde_test_repr[3:]) + pde_test_repr = EHField( + vector_from_device(queue, eval_repr_at(calc_patch_tgt))) maxwell_residuals = [ - calc_patch.norm(x, np.inf) / calc_patch.norm(pde_test_e, np.inf) + calc_patch.norm(x, np.inf) / calc_patch.norm(pde_test_repr.e, np.inf) for x in frequency_domain_maxwell( - calc_patch, pde_test_e, pde_test_h, case.k)] + calc_patch, pde_test_repr.e, pde_test_repr.h, case.k)] print("Maxwell residuals:", maxwell_residuals) eoc_rec_repr_maxwell.add_data_point(h_max, max(maxwell_residuals)) @@ -284,8 +284,8 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ ("j", jxyz), ("rho", rho), - ("Einc", inc_field_scat[:3]), - ("Hinc", inc_field_scat[3:]), + ("Einc", inc_field_scat.e), + ("Hinc", inc_field_scat.h), ("bdry_normals", bdry_normals), ]) @@ -310,15 +310,18 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): ]) raise - fplot_inc = eval_inc_field_at(fplot_tgt) + fplot_repr = EHField(vector_from_device(queue, fplot_repr)) + + fplot_inc = EHField( + vector_from_device(queue, eval_inc_field_at(fplot_tgt))) fplot.write_vtk_file( "potential-%s.vts" % resolution, [ - ("E", vector_from_device(queue, fplot_repr[:3])), - ("H", vector_from_device(queue, fplot_repr[3:])), - ("Einc", vector_from_device(queue, fplot_inc[:3])), - ("Hinc", vector_from_device(queue, fplot_inc[3:])), + ("E", fplot_repr.e), + ("H", fplot_repr.h), + ("Einc", fplot_inc.e), + ("Hinc", fplot_inc.h), ] ) @@ -326,22 +329,15 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): # {{{ error in E, H - obs_repr = eval_repr_at(obs_discr) - - obs_e = obs_repr[:3] - obs_h = obs_repr[3:] - - inc_obs_e = inc_field_obs[:3] - inc_obs_h = inc_field_obs[3:] + obs_repr = EHField(eval_repr_at(obs_discr)) def obs_norm(f): return norm(obs_discr, queue, f, p=np.inf) - # FIXME: Why "+"? - rel_err_e = (obs_norm(inc_obs_e + obs_e) - / obs_norm(inc_obs_e)) - rel_err_h = (obs_norm(inc_obs_h + obs_h) - / obs_norm(inc_obs_h)) + rel_err_e = (obs_norm(inc_field_obs.e - obs_repr.e) + / obs_norm(inc_field_obs.e)) + rel_err_h = (obs_norm(inc_field_obs.h - obs_repr.h) + / obs_norm(inc_field_obs.h)) # }}} -- GitLab From a0851a399343a297f5393d487ae351db65e874f1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 16 Sep 2017 19:59:49 -0500 Subject: [PATCH 36/53] Add test for tangential ONB functionality --- test/test_symbolic.py | 50 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 48 insertions(+), 2 deletions(-) diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 83f78d64..6894e15e 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -26,6 +26,10 @@ import pytest import numpy as np import pyopencl as cl +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +from pytential import bind, sym import logging logger = logging.getLogger(__name__) @@ -99,18 +103,21 @@ def get_unit_sphere_with_ref_mean_curvature(cl_ctx): # }}} +# {{{ test_mean_curvature + @pytest.mark.parametrize(("discr_name", "discr_and_ref_mean_curvature_getter"), [ ("unit_circle", get_ellipse_with_ref_mean_curvature), ("2-to-1 ellipse", partial(get_ellipse_with_ref_mean_curvature, aspect=2)), ("square", get_square_with_ref_mean_curvature), ("unit sphere", get_unit_sphere_with_ref_mean_curvature), ]) -def test_mean_curvature(ctx_getter, discr_name, discr_and_ref_mean_curvature_getter): +def test_mean_curvature(ctx_factory, discr_name, + discr_and_ref_mean_curvature_getter): if discr_name == "unit sphere": pytest.skip("not implemented in 3D yet") import pytential.symbolic.primitives as prim - ctx = ctx_getter() + ctx = ctx_factory() discr, ref_mean_curvature = discr_and_ref_mean_curvature_getter(ctx) @@ -122,6 +129,45 @@ def test_mean_curvature(ctx_getter, discr_name, discr_and_ref_mean_curvature_get assert np.allclose(mean_curvature.get(), ref_mean_curvature) +# }}} + + +# {{{ test_tangential_onb + +def test_tangential_onb(ctx_factory): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + + from meshmode.mesh.generation import generate_torus + mesh = generate_torus(5, 2, order=3) + + discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(3)) + + tob = sym.tangential_onb(mesh.ambient_dim) + nvecs = tob.shape[1] + + # make sure tangential_onb is mutually orthogonal and normalized + orth_check = bind(discr, sym.make_obj_array([ + np.dot(tob[:, i], tob[:, j]) - (1 if i == j else 0) + for i in range(nvecs) for j in range(nvecs)]) + )(queue) + + for i, orth_i in enumerate(orth_check): + assert (cl.clmath.fabs(orth_i) < 1e-13).get().all() + + # make sure tangential_onb is orthogonal to normal + orth_check = bind(discr, sym.make_obj_array([ + np.dot(tob[:, i], sym.normal(mesh.ambient_dim).as_vector()) + for i in range(nvecs)]) + )(queue) + + for i, orth_i in enumerate(orth_check): + assert (cl.clmath.fabs(orth_i) < 1e-13).get().all() + +# }}} + # You can test individual routines by typing # $ python test_tools.py 'test_routine()' -- GitLab From c7ff79ef8e21c1d278e915a42187af744ce0e1cc Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 12:36:55 -0500 Subject: [PATCH 37/53] Add test for 3D jump relations --- test/test_layer_pot.py | 147 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 146 insertions(+), 1 deletion(-) diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index d3c7e575..af642093 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -36,7 +36,7 @@ from meshmode.mesh.generation import ( # noqa ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle, make_curve_mesh, NArmedStarfish) from sumpy.visualization import FieldPlotter -from pytential import bind, sym +from pytential import bind, sym, norm import logging logger = logging.getLogger(__name__) @@ -386,6 +386,151 @@ def test_perf_data_gathering(ctx_getter, n_arms=5): # }}} +# {{{ test 3D jump relations + +@pytest.mark.parametrize("relation", ["sp", "nxcurls"]) +def test_3d_jump_relations(ctx_factory, relation, visualize=False): + #logging.basicConfig(level=logging.INFO) + + pytest.importorskip("pyfmmlib") + + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + + target_order = 4 + qbx_order = target_order + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + for nel_factor in [6, 8, 12]: + from meshmode.mesh.generation import generate_torus + mesh = generate_torus( + 5, 2, order=target_order, + n_outer=2*nel_factor, n_inner=nel_factor) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + pre_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(3)) + + from pytential.qbx import QBXLayerPotentialSource + qbx, _ = QBXLayerPotentialSource( + pre_discr, fine_order=4*target_order, + qbx_order=qbx_order, + fmm_order=qbx_order + 5, + fmm_backend="fmmlib" + ).with_refinement() + + from sumpy.kernel import LaplaceKernel + knl = LaplaceKernel(3) + + def nxcurlS(qbx_forced_limit): + + return sym.n_cross(sym.curl(sym.S( + knl, + sym.cse(sym.tangential_to_xyz(density_sym), "jxyz"), + qbx_forced_limit=qbx_forced_limit))) + + x, y, z = qbx.density_discr.nodes().with_queue(queue) + m = cl.clmath + + if relation == "nxcurls": + density_sym = sym.make_sym_vector("density", 2) + + jump_identity_sym = ( + nxcurlS(+1) + - (nxcurlS("avg") + 0.5*sym.tangential_to_xyz(density_sym))) + + # The tangential coordinate system is element-local, so we can't just + # conjure up some globally smooth functions, interpret their values + # in the tangential coordinate system, and be done. Instead, generate + # an XYZ function and project it. + density = bind( + qbx, + sym.xyz_to_tangential(sym.make_sym_vector("jxyz", 3)))( + queue, + jxyz=sym.make_obj_array([ + m.cos(0.5*x) * m.cos(0.5*y) * m.cos(0.5*z), + m.sin(0.5*x) * m.cos(0.5*y) * m.sin(0.5*z), + m.sin(0.5*x) * m.cos(0.5*y) * m.cos(0.5*z), + ])) + + elif relation == "sp": + + density = m.cos(2*x) * m.cos(2*y) * m.cos(z) + density_sym = sym.var("density") + + jump_identity_sym = ( + sym.Sp(knl, density_sym, qbx_forced_limit=+1) + - (sym.Sp(knl, density_sym, qbx_forced_limit="avg") + - 0.5*density_sym)) + + else: + raise ValueError("unexpected value of 'relation': %s" % relation) + + bound_jump_identity = bind(qbx, jump_identity_sym) + jump_identity = bound_jump_identity(queue, density=density) + + err = ( + norm(qbx, queue, jump_identity, np.inf) + / + norm(qbx, queue, density, np.inf)) + print("ERROR", qbx.h_max, err) + + eoc_rec.add_data_point(qbx.h_max, err) + + # {{{ visualization + + if visualize and relation == "nxcurls": + nxcurlS_ext = bind(qbx, nxcurlS(+1))(queue, density=density) + nxcurlS_avg = bind(qbx, nxcurlS("avg"))(queue, density=density) + jtxyz = bind(qbx, sym.tangential_to_xyz(density_sym))( + queue, density=density) + + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, qbx.density_discr, target_order+3) + + bdry_normals = bind(qbx, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("source-%s.vtu" % nel_factor, [ + ("jt", jtxyz), + ("nxcurlS_ext", nxcurlS_ext), + ("nxcurlS_avg", nxcurlS_avg), + ("bdry_normals", bdry_normals), + ]) + + if visualize and relation == "sp": + sp_ext = bind(qbx, sym.Sp(knl, density_sym, qbx_forced_limit=+1))( + queue, density=density) + sp_avg = bind(qbx, sym.Sp(knl, density_sym, qbx_forced_limit="avg"))( + queue, density=density) + + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, qbx.density_discr, target_order+3) + + bdry_normals = bind(qbx, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("source-%s.vtu" % nel_factor, [ + ("density", density), + ("sp_ext", sp_ext), + ("sp_avg", sp_avg), + ("bdry_normals", bdry_normals), + ]) + + # }}} + + print(eoc_rec) + + assert eoc_rec.order_estimate() >= qbx_order - 1.5 + +# }}} + + # You can test individual routines by typing # $ python test_layer_pot.py 'test_routine()' -- GitLab From 470dd36f256be5e0cf8906df94d43a60e7823eef Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 12:37:38 -0500 Subject: [PATCH 38/53] Fix CSE call in sym.normal() to make sure normals get a CSE prefix of 'normal' --- pytential/symbolic/primitives.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index d00aef68..ea70e8b5 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -382,7 +382,8 @@ def normal(ambient_dim, dim=None, where=None): return cse( # Dorst Section 3.7.2 pder << pder.I.inv(), - cse_scope.DISCRETIZATION) + "normal", + scope=cse_scope.DISCRETIZATION) def mean_curvature(ambient_dim, dim=None, where=None): -- GitLab From 634c6a24c78e1d8c025090fe25f78945a5b803e7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 15:08:46 -0500 Subject: [PATCH 39/53] MFIE fully tested and working --- pytential/symbolic/pde/maxwell/__init__.py | 43 +++++++++++----------- test/test_maxwell.py | 43 ++++++++++++++++------ 2 files changed, 54 insertions(+), 32 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 030c4185..84175d89 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -35,13 +35,24 @@ cse = sym.cse # {{{ Charge-Current MFIE class PECAugmentedMFIEOperator: - r"""Magnetic Field Integral Equation operator, - under the assumption of no surface charges. + r"""Magnetic Field Integral Equation operator with PEC boundary + conditions, under the assumption of no surface charges. see :file:`contrib/notes/mfie.tm` + The arguments *loc* below decide between the exterior and the + interior MFIE. The "exterior" (loc=1) MFIE enforces a zero field + on the interior of the integration surface, whereas the "interior" MFIE + (loc=-1) enforces a zero field on the exterior. + Uses the sign convention :math:`\exp(-1 \omega t)` for the sinusoidal time dependency. + + .. automethod:: j_operator + .. automethod:: j_rhs + .. automethod:: rho_operator + .. automethod:: rho_rhs + .. automethod:: scattered_volume_field """ def __init__(self, k=sym.var("k")): @@ -52,7 +63,7 @@ class PECAugmentedMFIEOperator: def j_operator(self, loc, Jt): Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") return xyz_to_tangential( - (loc*0.5)*Jxyz - sym.n_cross( + loc*0.5*Jxyz - sym.n_cross( sym.curl(sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit="avg")))) @@ -60,32 +71,22 @@ class PECAugmentedMFIEOperator: return xyz_to_tangential(sym.n_cross(Hinc_xyz)) def rho_operator(self, loc, rho): - return (loc*0.5)*rho+sym.Sp(self.kernel, rho, k=self.k) + 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 (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 = sym.S(self.kernel, Jxyz, k=self.k) - phi = sym.S(self.kernel, rho, k=self.k) - - # use - n x n x v = v_tangential - - E_scat = 1j*self.k*A - sym.grad(3, phi) + 0.5*loc*rho - H_scat = sym.curl(sym.S(self.kernel, Jxyz, k=self.k)) + (loc*0.5)*Jxyz - - return sym.join_fields(E_scat, H_scat) + + 1j*self.k*sym.n_dot(sym.S( + self.kernel, Jxyz, k=self.k, + # continuous--qbx_forced_limit doesn't really matter + qbx_forced_limit="avg"))) - def scattered_volume_field(self, Jt, rho): + def scattered_volume_field(self, Jt, rho, qbx_forced_limit=None): Jxyz = sym.cse(sym.tangential_to_xyz(Jt), "Jxyz") - A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=None) - phi = sym.S(self.kernel, rho, k=self.k, qbx_forced_limit=None) + A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=qbx_forced_limit) + phi = sym.S(self.kernel, rho, k=self.k, qbx_forced_limit=qbx_forced_limit) E_scat = 1j*self.k*A - sym.grad(3, phi) H_scat = sym.curl(A) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index 26ec2e70..587b368f 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -70,7 +70,6 @@ class SphereTestCase(TestCase): from meshmode.mesh.generation import generate_icosphere if self.is_interior: - return generate_icosphere(5, target_order) else: return generate_icosphere(0.5, target_order) @@ -128,7 +127,11 @@ class EHField(object): tc_int, tc_ext ]) -def test_mfie_from_source(ctx_getter, case, visualize=False): +def test_pec_mfie_extinction(ctx_getter, case, visualize=False): + """For (say) is_interior=False (the 'exterior' MFIE), this test verifies + extinction of the combined (incoming + scattered) field on the interior + of the scatterer. + """ logging.basicConfig(level=logging.INFO) cl_ctx = ctx_getter() @@ -176,14 +179,11 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): # }}} loc_sign = -1 if case.is_interior else +1 - # xyz_mfie_bdry_vol_op = mfie.boundary_field(0, Jxyz) - # xyz_mfie_vol_op = mfie.volume_field(Jxyz) - - # n_cross_op = n_cross(make_vector_field("vec", 3)) from pytools.convergence import EOCRecorder eoc_rec_repr_maxwell = EOCRecorder() + eoc_pec_bc = EOCRecorder() eoc_rec_e = EOCRecorder() eoc_rec_h = EOCRecorder() @@ -272,6 +272,29 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): # }}} + # {{{ check PEC BC on total field + + bc_repr = EHField(mfie.scattered_volume_field( + jt_sym, rho_sym, qbx_forced_limit=loc_sign)) + pec_bc_e = sym.n_cross(bc_repr.e + inc_xyz_sym.e) + pec_bc_h = sym.normal(3).as_vector().dot(bc_repr.h + inc_xyz_sym.h) + + eh_bc_values = bind(qbx, sym.join_fields(pec_bc_e, pec_bc_h))( + queue, jt=jt, rho=rho, inc_fld=inc_field_scat.field, + **knl_kwargs) + + def scat_norm(f): + return norm(qbx, queue, f, p=np.inf) + + e_bc_residual = scat_norm(eh_bc_values[:3]) / scat_norm(inc_field_scat.e) + h_bc_residual = scat_norm(eh_bc_values[3]) / scat_norm(inc_field_scat.h) + + print("E/H PEC BC residuals:", h_max, e_bc_residual, h_bc_residual) + + eoc_pec_bc.add_data_point(h_max, max(e_bc_residual, h_bc_residual)) + + # }}} + # {{{ visualization if visualize: @@ -334,9 +357,9 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): def obs_norm(f): return norm(obs_discr, queue, f, p=np.inf) - rel_err_e = (obs_norm(inc_field_obs.e - obs_repr.e) + rel_err_e = (obs_norm(inc_field_obs.e + obs_repr.e) / obs_norm(inc_field_obs.e)) - rel_err_h = (obs_norm(inc_field_obs.h - obs_repr.h) + rel_err_h = (obs_norm(inc_field_obs.h + obs_repr.h) / obs_norm(inc_field_obs.h)) # }}} @@ -346,9 +369,6 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): eoc_rec_h.add_data_point(h_max, rel_err_h) eoc_rec_e.add_data_point(h_max, rel_err_e) - # TODO: Check that total field verifies BC - # TODO: Check for extinction on the interior - print("--------------------------------------------------------") print("is_interior=%s" % case.is_interior) print("--------------------------------------------------------") @@ -356,6 +376,7 @@ def test_mfie_from_source(ctx_getter, case, visualize=False): good = True for which_eoc, eoc_rec, order_tol in [ ("maxwell", eoc_rec_repr_maxwell, 1.5), + ("PEC BC", eoc_pec_bc, 1), ("H", eoc_rec_h, 1), ("E", eoc_rec_e, 1)]: print(which_eoc) -- GitLab From bb18fe150c98f1f6815f983e31727f1e87d246da Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 15:29:37 -0500 Subject: [PATCH 40/53] Fix/update MFIE derivation notes --- contrib/notes/mfie.tm | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/contrib/notes/mfie.tm b/contrib/notes/mfie.tm index b2f2e648..110152dc 100644 --- a/contrib/notes/mfie.tm +++ b/contrib/notes/mfie.tm @@ -48,12 +48,14 @@ The interior MFIE is derived similarly as: <\equation*> - n\H=-|2>-\SJ|)>|)>. + n\H=-|2>-\SJ|)>|)> + (Note the sign flip originating from the boundary condition.) + > postprocessor derivation> - We'll start with the boundary condition + Start with the boundary condition <\equation*> n\-E|)>=\. @@ -65,23 +67,23 @@ n\+E|)>=\. - Now use the representation =-i*k*A-\\=i*k*SJ-\S\> + Now use the representation =i*k*A-\\=i*k*SJ-\S\> and obtain <\equation*> - n\E+J-\S\|)>|)>=\. + n\E+J-\S\|)>|)>=\. Carrying out the limit, we obtain: <\equation*> - n\E-n\J|)>-S\+\=\. + n\E+n\J|)>-S\+\=\. Rearrange: <\equation*> - n\E-n\J|)>=\+S\. + n\E+n\J|)>=\+S\. -- GitLab From 0bece3c56a6cff8c44214aa4efef8c9be6403d67 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 17:50:22 -0500 Subject: [PATCH 41/53] Add MFIE rounded-cube test case (not run by default) --- test/test_maxwell.py | 47 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index 587b368f..35b30887 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -89,11 +89,58 @@ class SphereTestCase(TestCase): cl.array.to_device(queue, sources)) +class RoundedCubeTestCase(TestCase): + target_order = 8 + gmres_tol = 1e-10 + + def get_mesh(self, resolution, target_order): + from meshmode.mesh.io import generate_gmsh, FileSource + mesh = generate_gmsh( + FileSource("rounded-cube.step"), 2, order=3, + other_options=[ + "-string", + "Mesh.CharacteristicLengthMax = %g;" % resolution]) + + from meshmode.mesh.processing import perform_flips, affine_map + mesh = affine_map(mesh, b=np.array([-0.5, -0.5, -0.5])) + mesh = affine_map(mesh, A=np.eye(3)*2) + + # now centered at origin and extends to -1,1 + + # Flip elements--gmsh generates inside-out geometry. + return perform_flips(mesh, np.ones(mesh.nelements)) + + def get_observation_mesh(self, target_order): + from meshmode.mesh.generation import generate_icosphere + + if self.is_interior: + return generate_icosphere(5, target_order) + else: + return generate_icosphere(0.5, target_order) + + def get_source(self, queue): + if self.is_interior: + source_ctr = np.array([[0.35, 0.1, 0.15]]).T + else: + source_ctr = np.array([[5, 0.1, 0.15]]).T + + source_rad = 0.3 + + sources = source_ctr + source_rad*2*(np.random.rand(3, 10)-0.5) + from pytential.source import PointPotentialSource + return PointPotentialSource( + queue.context, + cl.array.to_device(queue, sources)) + + tc_int = SphereTestCase(k=1.2, is_interior=True, resolutions=[0, 1], qbx_order=3, fmm_order=10) tc_ext = SphereTestCase(k=1.2, is_interior=False, resolutions=[0, 1], qbx_order=3, fmm_order=10) +tc_rc_ext = RoundedCubeTestCase(k=1.2, is_interior=False, resolutions=[0.3], + qbx_order=3, fmm_order=10) + def get_sym_maxwell_source(kernel, jxyz, k): # This ensures div A = 0, which is simply a consequence of div curl S=0. -- GitLab From c85572975320bac6ebe2d256672cdb91e719ba01 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 17:58:45 -0500 Subject: [PATCH 42/53] Bump down EOC bars for interior MFIE --- test/test_maxwell.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index 35b30887..c22457f7 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -423,9 +423,9 @@ def test_pec_mfie_extinction(ctx_getter, case, visualize=False): good = True for which_eoc, eoc_rec, order_tol in [ ("maxwell", eoc_rec_repr_maxwell, 1.5), - ("PEC BC", eoc_pec_bc, 1), - ("H", eoc_rec_h, 1), - ("E", eoc_rec_e, 1)]: + ("PEC BC", eoc_pec_bc, 1.5), + ("H", eoc_rec_h, 1.5), + ("E", eoc_rec_e, 1.5)]: print(which_eoc) print(eoc_rec.pretty_print()) -- GitLab From 58f114ed0ffc5521d31768246525ed970c08223f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 18:36:19 -0500 Subject: [PATCH 43/53] Make sym functions (e.g. sqrt) work on scalars (Closes #54 on gitlab) --- pytential/symbolic/execution.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 53826668..869978f7 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -146,8 +146,14 @@ class EvaluationMapper(EvaluationMapperBase): if isinstance(expr.function, EvalMapperFunction): return getattr(self, "apply_"+expr.function.name)(expr.parameters) elif isinstance(expr.function, CLMathFunction): - return getattr(cl.clmath, expr.function.name)( - *(self.rec(arg) for arg in expr.parameters), queue=self.queue) + args = [self.rec(arg) for arg in expr.parameters] + from numbers import Number + if all(isinstance(arg, Number) for arg in args): + return getattr(np, expr.function.name)(*args) + else: + return getattr(cl.clmath, expr.function.name)( + *args, queue=self.queue) + else: return EvaluationMapperBase.map_call(self, expr) -- GitLab From 72fa0176dc9008c03d76f3c98d41da6f70424f89 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 18:37:53 -0500 Subject: [PATCH 44/53] Allow bind(PointsTarget, ...) --- pytential/symbolic/execution.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 869978f7..3e96b00c 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -313,13 +313,14 @@ def prepare_places(places): from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET from meshmode.discretization import Discretization from pytential.source import LayerPotentialSourceBase + from pytential.target import TargetBase if isinstance(places, LayerPotentialSourceBase): places = { DEFAULT_SOURCE: places, DEFAULT_TARGET: places.density_discr, } - elif isinstance(places, Discretization): + elif isinstance(places, (Discretization, TargetBase)): places = { DEFAULT_TARGET: places, } -- GitLab From edb70978696f706001b274018fdb270c531a3e24 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 18:38:39 -0500 Subject: [PATCH 45/53] Primitives: Improve docs, add sym.div --- pytential/symbolic/primitives.py | 76 +++++++++++++++++++++++++++----- 1 file changed, 64 insertions(+), 12 deletions(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index ea70e8b5..a52f384b 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -86,9 +86,30 @@ Functions .. data:: real .. data:: imag .. data:: conj -.. data:: sqrt .. data:: abs +.. data:: sqrt + +.. data:: sin +.. data:: cos +.. data:: tan + +.. data:: asin +.. data:: acos +.. data:: atan +.. data:: atan2 + +.. data:: sinh +.. data:: cosh +.. data:: tanh + +.. data:: asinh +.. data:: acosh +.. data:: atanh + +.. data:: exp +.. data:: log + Discretization properties ^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -106,14 +127,17 @@ Elementary numerics .. autoclass:: NumReferenceDerivative .. autoclass:: NodeSum +.. autoclass:: NodeMax .. autofunction:: integral .. autoclass:: Ones .. autofunction:: ones_vec .. autofunction:: area +.. autofunction:: mean .. autoclass:: IterativeInverse -Calculus -^^^^^^^^ +Calculus (based on Geometric Algebra) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + .. autoclass:: Derivative .. autofunction:: dd_axis .. autofunction:: d_dx @@ -121,12 +145,39 @@ Calculus .. autofunction:: d_dz .. autofunction:: grad_mv .. autofunction:: grad +.. autofunction:: laplace Layer potentials ^^^^^^^^^^^^^^^^ .. autoclass:: IntG .. autofunction:: int_g_dsource + +.. autofunction:: S +.. autofunction:: Sp +.. autofunction:: Spp +.. autofunction:: D +.. autofunction:: Dp +.. autofunction:: normal_derivative +.. autofunction:: tangential_derivative + +"Conventional" Vector Calculus +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. autofunction:: tangential_onb +.. autofunction:: xyz_to_tangential +.. autofunction:: tangential_to_xyz +.. autofunction:: project_to_tangential +.. autofunction:: cross +.. autofunction:: n_dot +.. autofunction:: n_cross +.. autofunction:: curl +.. autofunction:: pretty + +Pretty-printing expressions +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. autofunction:: pretty """ @@ -961,7 +1012,7 @@ def Dp(kernel, *args, **kwargs): # noqa # }}} -# {{{ geometric operations +# {{{ conventional vector calculus def tangential_onb(ambient_dim, dim=None, where=None): if dim is None: @@ -1015,24 +1066,25 @@ def n_dot(vec, where=None): return np.dot(nrm, vec) -def n_cross(vec, where=None): - nrm = normal(3, where).as_vector() +def cross(vec_a, vec_b): + assert len(vec_a) == len(vec_b) == 3 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] + levi_civita((i, j, k)) * vec_a[j] * vec_b[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) +def n_cross(vec, where=None): + return cross(normal(3, where).as_vector(), 3) + - return sum(dd_axis(iaxis, ambient_dim, S(kernel, arg[iaxis])) - for iaxis in range(ambient_dim)) +def div(vec): + ambient_dim = len(vec) + return sum(dd_axis(iaxis, ambient_dim, vec) for iaxis in range(ambient_dim)) def curl(vec): -- GitLab From bd72bee35e81cdb779b755d4e43abdb140788be5 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 18:39:09 -0500 Subject: [PATCH 46/53] Scalar PDE ops: doc improvements --- pytential/symbolic/pde/scalar.py | 40 ++++++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 10 deletions(-) diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py index 15f542b7..fc64b839 100644 --- a/pytential/symbolic/pde/scalar.py +++ b/pytential/symbolic/pde/scalar.py @@ -23,16 +23,9 @@ THE SOFTWARE. """ __doc__ = """ - .. autoclass:: L2WeightedPDEOperator .. autoclass:: DirichletOperator .. autoclass:: NeumannOperator - -2D Dielectric -^^^^^^^^^^^^^ - -.. autoclass:: DielectricSRep2DBoundaryOperator -.. autoclass:: DielectricSDRep2DBoundaryOperator """ @@ -76,9 +69,20 @@ class L2WeightedPDEOperator(object): # {{{ dirichlet class DirichletOperator(L2WeightedPDEOperator): - """When testing this as a potential matcher, note that it can only - access potentials that come from charge distributions having *no* net - charge. (This is true at least in 2D.) + """IE operator and field representation for solving Dirichlet boundary + value problems with scalar kernels (e.g. + :class:`sumpy.kernel.LaplaceKernel`, + :class:`sumpy.kernel.HelmholtzKernel`, :class:`sumpy.kernel.YukawaKernel`) + + ..note :: + + When testing this as a potential matcher, note that it can only + access potentials that come from charge distributions having *no* net + charge. (This is true at least in 2D.) + + .. automethod:: is_unique_only_up_to_constant + .. automethod:: representation + .. automethod:: operator """ def __init__(self, kernel, loc_sign, alpha=None, use_l2_weighting=False, @@ -171,6 +175,22 @@ class DirichletOperator(L2WeightedPDEOperator): # {{{ neumann class NeumannOperator(L2WeightedPDEOperator): + """IE operator and field representation for solving Dirichlet boundary + value problems with scalar kernels (e.g. + :class:`sumpy.kernel.LaplaceKernel`, + :class:`sumpy.kernel.HelmholtzKernel`, :class:`sumpy.kernel.YukawaKernel`) + + ..note :: + + When testing this as a potential matcher, note that it can only + access potentials that come from charge distributions having *no* net + charge. (This is true at least in 2D.) + + .. automethod:: is_unique_only_up_to_constant + .. automethod:: representation + .. automethod:: operator + """ + def __init__(self, kernel, loc_sign, alpha=None, use_improved_operator=True, laplace_kernel=0, use_l2_weighting=False, -- GitLab From 38f77c1644172ef2789a3126c409acb81280796a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 18:40:03 -0500 Subject: [PATCH 47/53] Move point source function to pde.maxwell, add plane wave source, fix docs --- pytential/symbolic/pde/maxwell/__init__.py | 82 +++++++++++++++++++++- test/test_maxwell.py | 44 ++++++------ 2 files changed, 103 insertions(+), 23 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 84175d89..bcef1ee4 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -31,10 +31,90 @@ tangential_to_xyz = sym.tangential_to_xyz xyz_to_tangential = sym.xyz_to_tangential cse = sym.cse +__doc__ = """ + +.. autofunction:: get_sym_maxwell_point_source +.. autofunction:: get_sym_maxwell_plane_wave +.. autoclass:: PECChargeCurrentMFIEOperator +""" + + +# {{{ point source + +def get_sym_maxwell_point_source(kernel, jxyz, k): + """Return a symbolic expression that, when bound to a + :class:`pytential.source.PointPotentialSource` will yield + a field satisfying Maxwell's equations. + + Uses the sign convention :math:`\exp(-1 \omega t)` for the time dependency. + + This will return an object of six entries, the first three of which + represent the electric, and the second three of which represent the + magnetic field. This satisfies the time-domain Maxwell's equations + as verified by :func:`sumpy.point_calculus.frequency_domain_maxwell`. + """ + # This ensures div A = 0, which is simply a consequence of div curl S=0. + # This means we use the Coulomb gauge to generate this field. + + A = sym.curl(sym.S(kernel, jxyz, k=k, qbx_forced_limit=None)) + + # https://en.wikipedia.org/w/index.php?title=Maxwell%27s_equations&oldid=798940325#Alternative_formulations + # (Vector calculus/Potentials/Any Gauge) + # assumed time dependence exp(-1j*omega*t) + return sym.join_fields( + 1j*k*A, + sym.curl(A)) + +# }}} + + +# {{{ plane wave + +def get_sym_maxwell_plane_wave(amplitude_vec, v, omega, epsilon=1, mu=1, where=None): + """Return a symbolic expression that, when bound to a + :class:`pytential.source.PointPotentialSource` will yield + a field satisfying Maxwell's equations. + + :arg amplitude_vec: should be orthogonal to *v*. If it is not, + it will be orthogonalized. + :arg v: a three-vector representing the phase velocity of the wave + (may be an object array of variables or a vector of concrete numbers) + While *v* may mathematically be complex-valued, this function + is for now only tested for real values. + :arg omega: Accepts the "Helmholtz k" to be compatible with other parts + of this module. + + Uses the sign convention :math:`\exp(-1 \omega t)` for the time dependency. + + This will return an object of six entries, the first three of which + represent the electric, and the second three of which represent the + magnetic field. This satisfies the time-domain Maxwell's equations + as verified by :func:`sumpy.point_calculus.frequency_domain_maxwell`. + """ + + # See section 7.1 of Jackson, third ed. for derivation. + + # NOTE: for complex, need to ensure real(n).dot(imag(n)) = 0 (7.15) + + x = sym.nodes(3, where).as_vector() + + v_mag_squared = sym.cse(np.dot(v, v), "v_mag_squared") + n = v/sym.sqrt(v_mag_squared) + + amplitude_vec = amplitude_vec - np.dot(amplitude_vec, n)*n + + c_inv = np.sqrt(mu*epsilon) + + e = amplitude_vec * sym.exp(1j*np.dot(n*omega, x)) + + return sym.join_fields(e, c_inv * sym.cross(n, e)) + +# }}} + # {{{ Charge-Current MFIE -class PECAugmentedMFIEOperator: +class PECChargeCurrentMFIEOperator: r"""Magnetic Field Integral Equation operator with PEC boundary conditions, under the assumption of no surface charges. diff --git a/test/test_maxwell.py b/test/test_maxwell.py index c22457f7..504eb9cb 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -142,20 +142,6 @@ tc_rc_ext = RoundedCubeTestCase(k=1.2, is_interior=False, resolutions=[0.3], qbx_order=3, fmm_order=10) -def get_sym_maxwell_source(kernel, jxyz, k): - # This ensures div A = 0, which is simply a consequence of div curl S=0. - # This means we use the Coulomb gauge to generate this field. - - A = sym.curl(sym.S(kernel, jxyz, k=k, qbx_forced_limit=None)) - - # https://en.wikipedia.org/w/index.php?title=Maxwell%27s_equations&oldid=798940325#Alternative_formulations - # (Vector calculus/Potentials/Any Gauge) - # assumed time dependence exp(-1j*omega*t) - return sym.join_fields( - 1j*k*A, - sym.curl(A)) - - class EHField(object): def __init__(self, eh_field): assert len(eh_field) == 6 @@ -196,22 +182,36 @@ def test_pec_mfie_extinction(ctx_getter, case, visualize=False): jt_sym = sym.make_sym_vector("jt", 2) rho_sym = sym.var("rho") - from pytential.symbolic.pde.maxwell import PECAugmentedMFIEOperator - mfie = PECAugmentedMFIEOperator() + from pytential.symbolic.pde.maxwell import ( + PECChargeCurrentMFIEOperator, + get_sym_maxwell_point_source, + get_sym_maxwell_plane_wave) + mfie = PECChargeCurrentMFIEOperator() test_source = case.get_source(queue) calc_patch = CalculusPatch(np.array([-3, 0, 0]), h=0.01) - calc_patch_tgt = PointsTarget(calc_patch.points) + calc_patch_tgt = PointsTarget(cl.array.to_device(queue, calc_patch.points)) rng = cl.clrandom.PhiloxGenerator(cl_ctx, seed=12) src_j = rng.normal(queue, (3, test_source.nnodes), dtype=np.float64) def eval_inc_field_at(tgt): - return bind( - (test_source, tgt), - get_sym_maxwell_source(mfie.kernel, j_sym, mfie.k) - )(queue, j=src_j, k=case.k) + if 0: + # plane wave + return bind( + tgt, + get_sym_maxwell_plane_wave( + amplitude_vec=np.array([1, 1, 1]), + v=np.array([1, 0, 0]), + omega=case.k) + )(queue) + else: + # point source + return bind( + (test_source, tgt), + get_sym_maxwell_point_source(mfie.kernel, j_sym, mfie.k) + )(queue, j=src_j, k=case.k) pde_test_inc = EHField( vector_from_device(queue, eval_inc_field_at(calc_patch_tgt))) @@ -369,7 +369,7 @@ def test_pec_mfie_extinction(ctx_getter, case, visualize=False): qbx_tgt_tol = qbx.copy(target_association_tolerance=0.2) - fplot_tgt = PointsTarget(fplot.points) + fplot_tgt = PointsTarget(cl.array.to_device(queue, fplot.points)) try: fplot_repr = eval_repr_at(fplot_tgt, source=qbx_tgt_tol) except QBXTargetAssociationFailedException as e: -- GitLab From 4f208c9e7d35698142d895469473e8873d332549 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 18:40:22 -0500 Subject: [PATCH 48/53] Maxwell module: minor doc improvements --- pytential/symbolic/pde/maxwell/__init__.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index bcef1ee4..489f97fe 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -118,14 +118,15 @@ class PECChargeCurrentMFIEOperator: r"""Magnetic Field Integral Equation operator with PEC boundary conditions, under the assumption of no surface charges. - see :file:`contrib/notes/mfie.tm` + See :file:`contrib/notes/mfie.tm` in the repository for a derivation. The arguments *loc* below decide between the exterior and the interior MFIE. The "exterior" (loc=1) MFIE enforces a zero field on the interior of the integration surface, whereas the "interior" MFIE (loc=-1) enforces a zero field on the exterior. - Uses the sign convention :math:`\exp(-1 \omega t)` + Uses the sign convention :math:`\exp(-1 \omega t)` for the time dependency. + for the sinusoidal time dependency. .. automethod:: j_operator @@ -163,6 +164,12 @@ class PECChargeCurrentMFIEOperator: qbx_forced_limit="avg"))) def scattered_volume_field(self, Jt, rho, qbx_forced_limit=None): + """ + This will return an object of six entries, the first three of which + represent the electric, and the second three of which represent the + magnetic field. This satisfies the time-domain Maxwell's equations + as verified by :func:`sumpy.point_calculus.frequency_domain_maxwell`. + """ Jxyz = sym.cse(sym.tangential_to_xyz(Jt), "Jxyz") A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=qbx_forced_limit) -- GitLab From d52ec52b8fa06c9e83898178de69df4461616ea7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 17 Sep 2017 23:49:22 -0500 Subject: [PATCH 49/53] Fix implementation of sym.n_cross() in terms of sym.cross() --- pytential/symbolic/primitives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index a52f384b..e47ce396 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -1079,7 +1079,7 @@ def cross(vec_a, vec_b): def n_cross(vec, where=None): - return cross(normal(3, where).as_vector(), 3) + return cross(normal(3, where).as_vector(), vec) def div(vec): -- GitLab From 7bfacee11eb14acd5cafba614e7b1d5822eba6b2 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 18 Sep 2017 10:13:55 -0500 Subject: [PATCH 50/53] Restrict pytential tests to run on CIs tagged large-node --- .gitlab-ci.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 956ebb97..334a41d0 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -8,6 +8,7 @@ Python 3.5 POCL: tags: - python3.5 - pocl + - large-node except: - tags @@ -21,6 +22,7 @@ Python 3.6 POCL: tags: - python3.6 - pocl + - large-node except: - tags @@ -33,6 +35,7 @@ Python 3.5 Conda: - ". ./build-and-test-py-project-within-miniconda.sh" tags: - linux + - large-node except: - tags @@ -46,6 +49,7 @@ Python 2.7 POCL: tags: - python2.7 - pocl + - large-node except: - tags -- GitLab From b4c928ed79833a8c22ed085a12840db1d04e4a00 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 18 Sep 2017 10:14:14 -0500 Subject: [PATCH 51/53] Disable interior MFIE test to save testing time --- test/test_maxwell.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index 504eb9cb..a96f3d71 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -157,8 +157,8 @@ class EHField(object): @pytest.mark.parametrize("case", [ - tc_int, - tc_ext + #tc_int, + tc_ext, ]) def test_pec_mfie_extinction(ctx_getter, case, visualize=False): """For (say) is_interior=False (the 'exterior' MFIE), this test verifies -- GitLab From 7799095666bbe05f37685438ab9bc816bf203c55 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 18 Sep 2017 11:51:04 -0500 Subject: [PATCH 52/53] MFIE test: skimp on FMM order to gain execution speed --- test/test_maxwell.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_maxwell.py b/test/test_maxwell.py index a96f3d71..4228f4df 100644 --- a/test/test_maxwell.py +++ b/test/test_maxwell.py @@ -134,9 +134,9 @@ class RoundedCubeTestCase(TestCase): tc_int = SphereTestCase(k=1.2, is_interior=True, resolutions=[0, 1], - qbx_order=3, fmm_order=10) + qbx_order=3, fmm_order=5) tc_ext = SphereTestCase(k=1.2, is_interior=False, resolutions=[0, 1], - qbx_order=3, fmm_order=10) + qbx_order=3, fmm_order=5) tc_rc_ext = RoundedCubeTestCase(k=1.2, is_interior=False, resolutions=[0.3], qbx_order=3, fmm_order=10) -- GitLab From f4399df1267fa2922e2aa0bd5415d047327f617a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 18 Sep 2017 13:39:46 -0500 Subject: [PATCH 53/53] Add explanatory note to (unfinished) examples/maxwell.py --- examples/maxwell.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/maxwell.py b/examples/maxwell.py index a88d1bdf..f1c57e95 100644 --- a/examples/maxwell.py +++ b/examples/maxwell.py @@ -1,3 +1,6 @@ +# This is an untested/non-working work in progress. For a PEC Maxwell solver, +# see test/test_maxwell.py. + from __future__ import division import numpy as np import pyopencl as cl -- GitLab