diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 956ebb97b18040c87be62149560416b1ab9a5d41..334a41d0fb0b56e6582d71e9242c09d9729fdb0e 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 diff --git a/contrib/notes/.gitignore b/contrib/notes/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..55acb67b8aa39c019f7b7ec86ab21780bd927e02 --- /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 0000000000000000000000000000000000000000..110152dc94be3f01be6016ae821c96621463ac54 --- /dev/null +++ b/contrib/notes/mfie.tm @@ -0,0 +1,117 @@ + + + + +<\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|)>|)> + + + (Note the sign flip originating from the boundary condition.) + + > postprocessor derivation> + + 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||1Augmented + MFIE derivation> |.>>>>|> + + + |1.1MFIE derivation + |.>>>>|> + > + + |1.2|\> + postprocessor derivation |.>>>>|> + > + + + \ No newline at end of file diff --git a/examples/maxwell.py b/examples/maxwell.py index ed8c6106fcde56e7919344fb44b9500f7013aed0..f1c57e95718fdff437b61cd4b975f9bda668908a 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 @@ -18,7 +21,6 @@ queue = cl.CommandQueue(cl_ctx) target_order = 5 qbx_order = 3 -nelements = 60 mode_nr = 0 h = 1 @@ -27,6 +29,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, @@ -56,31 +59,176 @@ 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=0.4, + epss=[1.4, 1.0], + mus=[1.2, 1.0], + ) 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) - # {{{ make a density + if 1: + expr = sym_repr + print(sym.pretty(expr)) - nodes = density_discr.nodes().with_queue(queue) + print("#"*80) + from pytential.target import PointsTarget - angle = cl.clmath.atan2(nodes[1], nodes[0]) + tgt_points=np.zeros((3,1)) + tgt_points[0,0] = 100 + tgt_points[1,0] = -200 + tgt_points[2,0] = 300 - 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 + bound_op = bind((qbx, PointsTarget(tgt_points)), expr) + print(bound_op.code) - sigma = sigma.astype(np.complex128) + if 1: - jt = sym.make_obj_array([sigma, sigma]) - rho = sigma + 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 + 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 + + 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 # }}} @@ -102,6 +250,9 @@ def main(): qbx_tgt_tol = qbx.copy(target_association_tolerance=0.1) from pytential.target import PointsTarget from pytential.qbx import QBXTargetAssociationFailedException + + rho_sym = sym.var("rho") + try: fld_in_vol = bind( (qbx_tgt_tol, PointsTarget(fplot.points)), diff --git a/examples/maxwell_sphere.py b/examples/maxwell_sphere.py new file mode 100644 index 0000000000000000000000000000000000000000..d0cf8b3dd14015f7fca2e883e4acfb64d24740f0 --- /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/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index f146f2b694f47c0aeb6ca3ffbe0ec5d3f9640bc0..ebadacc80fada728d2782cc92ad1bb95daaf359c 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) diff --git a/pytential/source.py b/pytential/source.py index f6cabf2087859686ce1a6975c51c6c26d3e8d267..7e4f81fec77eb38ca24bb89e5723e5fd4978a122 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -54,9 +54,11 @@ class PotentialSource(object): class PointPotentialSource(PotentialSource): """ - ... attributes:: points + ... 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 { diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 53826668d7d9d87599659a464d3dc12c29fbaab9..3e96b00c8ff6d951f346b13e4b01db48cb6a3c3a 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) @@ -307,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, } diff --git a/pytential/symbolic/old_diffop_primitives.py b/pytential/symbolic/old_diffop_primitives.py index 23028e14f635a6eb44b038742637578c00cc6d1c..66f85d8d5a1cc5b1e7954fc27c22393ca0c93773 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 diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 9096998a9663c5fe43150505595c413d6d0c287a..489f97fe8f34325a2c5a2bcc7f045bb596cc3775 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -22,435 +22,265 @@ 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 -# {{{ MFIE +__doc__ = """ -class PECAugmentedMFIEOperator: - """Magnetic Field Integral Equation operator, - under the assumption of no surface charges. - - see notes/mfie.tm - """ - - def __init__(self, k=sym.var("k")): - from sumpy.kernel import HelmholtzKernel - self.kernel = HelmholtzKernel(3) - self.k = k - - 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)) - - def j_rhs(self, Hinc_xyz): - return xyz_to_tangential(n_cross(Hinc_xyz)) - - def rho_operator(self, loc, rho): - return (loc*0.5)*rho+Sp(self.k, rho) - - 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)) - - 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) +.. autofunction:: get_sym_maxwell_point_source +.. autofunction:: get_sym_maxwell_plane_wave +.. autoclass:: PECChargeCurrentMFIEOperator +""" - # 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 +# {{{ point source - from pytools.obj_array import join_fields - return join_fields(E_scat, H_scat) +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. - def scattered_volume_field(self, Jt, rho): - Jxyz = sym.cse(sym.tangential_to_xyz(Jt), "Jxyz") + Uses the sign convention :math:`\exp(-1 \omega t)` for the time dependency. - A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=None) - #grad_phi = grad_S(self.kernel, rho, 3, k=self.k) + 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. - E_scat = 1j*self.k*A# - grad_phi - #H_scat = curl_S_volume(k, Jxyz) + A = sym.curl(sym.S(kernel, jxyz, k=k, qbx_forced_limit=None)) - from pytools.obj_array import join_fields - return E_scat #join_fields(E_scat, H_scat) + # 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)) # }}} -# {{{ generalized Debye +# {{{ plane wave -def _debye_S0_surf_div(nxnxE): - # Note integration by parts: S0 div_\Gamma (X) -> - \int \grad g0 \cdot +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. - nxnxE_t = cse(xyz_to_tangential(nxnxE), "nxnxE_t") + :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. - return - sum([ - IntGdSource(0, nxnxE_t[i], - ds_direction=make_tangent(i, 3)) - for i in range(3-1)]) + Uses the sign convention :math:`\exp(-1 \omega t)` for the time dependency. - -class DebyeOperatorBase(object): + 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 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] - # }}} + # See section 7.1 of Jackson, third ed. for derivation. - def m(self, j): - return n_cross(j) + # NOTE: for complex, need to ensure real(n).dot(imag(n)) = 0 (7.15) - def boundary_field_base(self, r, q, j): - k = self.k + x = sym.nodes(3, where).as_vector() - surf_grad_phi = surf_grad_S(self.k, r) + v_mag_squared = sym.cse(np.dot(v, v), "v_mag_squared") + n = v/sym.sqrt(v_mag_squared) - dpsi_dn = Sp(self.k, q) - dpsi_dn -= self.loc_sign*1/2*q + amplitude_vec = amplitude_vec - np.dot(amplitude_vec, n)*n - m = cse(self.m(j), "m") + c_inv = np.sqrt(mu*epsilon) - A = cse(S(k, j), "A") - Q = cse(S(k, m), "Q") + e = amplitude_vec * sym.exp(1j*np.dot(n*omega, x)) - nxnxE = ( - n_cross( - n_cross(1j * k * A) + return sym.join_fields(e, c_inv * sym.cross(n, e)) - - 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) +# {{{ Charge-Current MFIE - 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)) +class PECChargeCurrentMFIEOperator: + r"""Magnetic Field Integral Equation operator with PEC boundary + conditions, under the assumption of no surface charges. - def cluster_points(self): - return (-1/4, +1/2) + 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)` for the time dependency. + for the sinusoidal time dependency. + .. automethod:: j_operator + .. automethod:: j_rhs + .. automethod:: rho_operator + .. automethod:: rho_rhs + .. automethod:: scattered_volume_field + """ -class InvertingDebyeOperatorBase(DebyeOperatorBase): - def r_and_q(self, r, q): - return r, q + def __init__(self, k=sym.var("k")): + from sumpy.kernel import HelmholtzKernel + self.kernel = HelmholtzKernel(3) + self.k = k - 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 j_operator(self, loc, Jt): + Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") + return xyz_to_tangential( + loc*0.5*Jxyz - sym.n_cross( + sym.curl(sym.S(self.kernel, Jxyz, k=self.k, + qbx_forced_limit="avg")))) - 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 j_rhs(self, Hinc_xyz): + return xyz_to_tangential(sym.n_cross(Hinc_xyz)) - def integral_equation(self, r_tilde, q_tilde, hvf_coefficients): - fix = 0 + def rho_operator(self, loc, rho): + return loc*0.5*rho+sym.Sp(self.kernel, rho, k=self.k) - if self.invertible: - s_ones = cse(S(0,Ones()), "s_ones") + def rho_rhs(self, Jt, Einc_xyz): + Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") - def inv_rank_one_coeff(u): - return cse(Mean(u)) + return (sym.n_dot(Einc_xyz) + + 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, 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") - r_coeff = inv_rank_one_coeff(r_tilde) - q_coeff = inv_rank_one_coeff(q_tilde) + 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) - from pytools.obj_array import join_fields - factors = self.cluster_points() + E_scat = 1j*self.k*A - sym.grad(3, phi) + H_scat = sym.curl(A) - fix = join_fields( - factors[0]*s_ones*r_coeff, - factors[1]*Ones()*q_coeff, - ) + return sym.join_fields(E_scat, H_scat) - return DebyeOperatorBase.integral_equation( - self, r_tilde, q_tilde, hvf_coefficients, fix=fix) +# }}} +# {{{ 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) + self.omega = omega + self.mus = mus + self.epss = epss + self.ks = [ + sym.cse(omega*(eps*mu)**0.5, "k%d" % i) + for i, (eps, mu) in enumerate(zip(epss, mus))] -class InvertingDebyeOperator(InvertingDebyeOperatorBase): - "Debye operator based on r and q." + def make_unknown(self, name): + return sym.make_sym_vector(name, 6) - 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.) + unk_structure = namedtuple("MuellerUnknowns", ["jt", "rho_e", "mt", "rho_m"]) - 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")) + def split_unknown(self, unk): + return self.unk_structure( + jt=unk[:2], + rho_e=unk[2], + mt=unk[3:5], + rho_m=unk[5]) - 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") + def operator(self, unk): + u = self.split_unknown(unk) - return ( - 1j * self.k * ( - surf_grad_alpha_r - n_cross(surf_grad_alpha_q)) - + self.harmonic_vector_field_current(hvf_coefficients)) + Jxyz = cse(tangential_to_xyz(u.jt), "Jxyz") + Mxyz = cse(tangential_to_xyz(u.mt), "Mxyz") + omega = self.omega + mu0, mu1 = self.mus + eps0, eps1 = self.epss + k0, k1 = self.ks -class InvertingSLapSDebyeOperator(InvertingDebyeOperatorBase): - "Debye operator based on r and q." + S = partial(sym.S, self.kernel, qbx_forced_limit="avg") - 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.) + def curl_S(dens): + return sym.curl(sym.S(self.kernel, dens, qbx_forced_limit="avg")) - 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") + grad = partial(sym.grad, 3) - 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") + 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") - return ( - 1j * self.k * ( - surf_grad_alpha_r - n_cross(surf_grad_alpha_q)) - + self.harmonic_vector_field_current(hvf_coefficients)) + 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)) + # sign flip included + F4 = -sym.n_dot(mu1*H1-mu0*H0) + 0.5*(mu1+mu0)*u.rho_m -class NonInvertingDebyeOperator(DebyeOperatorBase): - "Debye operator based on r_tilde and q_tilde." + return sym.join_fields(F1, F2, F3, F4) - 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 rhs(self, Einc_xyz, Hinc_xyz): + mu1 = self.mus[1] + eps1 = self.epss[1] - 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 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)), + sym.n_dot(-mu1*Hinc_xyz)) - return (1j * self.k * ( - surf_grad_s2_r_tilde - n_cross(surf_grad_s2_q_tilde)) - + self.harmonic_vector_field_current(hvf_coefficients)) + 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 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) + # omega = self.omega + mu = self.mus[i] + eps = self.epss[i] + k = self.ks[i] - 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) + S = partial(sym.S, self.kernel, qbx_forced_limit=None, k=k) - def integral_equation(self, r_tilde, q_tilde, hvf_coefficients): - fix = 0 + def curl_S(dens): + return sym.curl(sym.S(self.kernel, dens, qbx_forced_limit=None, k=k)) - if self.invertible: - s_ones = cse(S(0, Ones()), "s_ones") + grad = partial(sym.grad, 3) - def inv_rank_one_coeff(u): - return cse(Mean(cse(S(0, cse(S(0, u)))))) + 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)) - 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) + return sym.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 0000000000000000000000000000000000000000..1b6ea8ee9dab508237973f1861ecb77405c495ab --- /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/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py index 15f542b764e767898c969839b2c3cd4dc973aede..fc64b839cd99b9ca29e88ac48181b6595cd20b04 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, diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 9869afe0f705a6c4b6113d68495b03c2af68e8c9..e47ce396baccee75a7646cd9470285cc00cc8af2 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 @@ -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 """ @@ -374,7 +425,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) @@ -382,7 +433,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): @@ -847,7 +899,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, @@ -960,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: @@ -977,8 +1029,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) # }}} @@ -1002,9 +1055,47 @@ 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 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)) * vec_a[j] * vec_b[k] + for j in range(3) for k in range(3)) + for i in range(3)]) + + +def n_cross(vec, where=None): + return cross(normal(3, where).as_vector(), vec) + + +def div(vec): + ambient_dim = len(vec) + return sum(dd_axis(iaxis, ambient_dim, vec) for iaxis in range(ambient_dim)) + + +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, vec[n]) + for m in range(3) for n in range(3)) + for l in range(3)]) # }}} diff --git a/setup.cfg b/setup.cfg index dee4ca77dc4e0f2cd8b9ebe3166f7954b392182d..42291e82433c3c7522f8e267e41321c4d19db099 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,7 +1,7 @@ [flake8] -ignore = E126,E127,E128,E123,E226,E241,E242,E265,E402,W503,D102,D103 +ignore = E126,E127,E128,E123,E226,E241,E242,E265,E402,W503,N803,N806,N802,D102,D103 max-line-length=85 exclude= - pytential/symbolic/pde/maxwell/__init__.py, pytential/symbolic/old_diffop_primitives.py, + pytential/symbolic/pde/maxwell/generalized_debye.py, diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index d3c7e5756e011b1b70484017f2519609a49ea868..af642093ebd40947038d17daf9b199658b06f105 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()' diff --git a/test/test_maxwell.py b/test/test_maxwell.py new file mode 100644 index 0000000000000000000000000000000000000000..4228f4dfc980750b4ece5f238462d152f66b8690 --- /dev/null +++ b/test/test_maxwell.py @@ -0,0 +1,450 @@ +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, norm + +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 +logger = logging.getLogger(__name__) + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + + +class TestCase: + fmm_backend = "fmmlib" + + 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): + 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 + return refine_uniformly( + generate_icosphere(2, target_order), + resolution) + + 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)) + + +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=5) +tc_ext = SphereTestCase(k=1.2, is_interior=False, resolutions=[0, 1], + 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) + + +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, + ]) +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() + queue = cl.CommandQueue(cl_ctx) + + pytest.importorskip("pyfmmlib") + + 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 ( + 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(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): + 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))) + + source_maxwell_resids = [ + 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_inc.e, pde_test_inc.h, case.k)] + print("Source Maxwell residuals:", source_maxwell_resids) + assert max(source_maxwell_resids) < 1e-6 + + # }}} + + loc_sign = -1 if case.is_interior else +1 + + from pytools.convergence import EOCRecorder + + eoc_rec_repr_maxwell = EOCRecorder() + eoc_pec_bc = EOCRecorder() + 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 = case.get_mesh(resolution, case.target_order) + observation_mesh = case.get_observation_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() + h_max = qbx.h_max + + scat_discr = qbx.density_discr + obs_discr = Discretization( + cl_ctx, observation_mesh, + InterpolatoryQuadratureSimplexGroupFactory(case.target_order)) + + inc_field_scat = EHField(eval_inc_field_at(scat_discr)) + inc_field_obs = EHField(eval_inc_field_at(obs_discr)) + + # {{{ system solve + + 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(inc_xyz_sym.h))( + queue, inc_fld=inc_field_scat.field, **knl_kwargs) + + 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) + + 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, 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), + rho_rhs, + tol=case.gmres_tol, + progress=True, + hard_failure=True) + + rho = gmres_result.solution + + # }}} + + 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, source=None): + if source is None: + source = qbx + + return bind((source, tgt), sym_repr)(queue, jt=jt, rho=rho, **knl_kwargs) + + 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_repr.e, np.inf) + for x in frequency_domain_maxwell( + 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)) + + # }}} + + # {{{ 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: + 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.e), + ("Hinc", inc_field_scat.h), + ("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.2) + + 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: + fplot.write_vtk_file( + "failed-targets.vts", + [ + ("failed_targets", e.failed_target_flags.get(queue)) + ]) + raise + + 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", fplot_repr.e), + ("H", fplot_repr.h), + ("Einc", fplot_inc.e), + ("Hinc", fplot_inc.h), + ] + ) + + # }}} + + # {{{ error in E, H + + obs_repr = EHField(eval_repr_at(obs_discr)) + + 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) + / obs_norm(inc_field_obs.e)) + rel_err_h = (obs_norm(inc_field_obs.h + obs_repr.h) + / obs_norm(inc_field_obs.h)) + + # }}} + + 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) + + print("--------------------------------------------------------") + print("is_interior=%s" % case.is_interior) + print("--------------------------------------------------------") + + good = True + for which_eoc, eoc_rec, order_tol in [ + ("maxwell", eoc_rec_repr_maxwell, 1.5), + ("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()) + + if len(eoc_rec.history) > 1: + if eoc_rec.order_estimate() < case.qbx_order - order_tol: + good = False + + assert good + + +# 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 diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py index ec8a4c94c569843228e451948bf7b8231d599e25..690541a5f9df173f9b2df688c038e9b44c31a83d 100644 --- a/test/test_scalar_int_eq.py +++ b/test/test_scalar_int_eq.py @@ -605,6 +605,35 @@ class EllipsoidIntEqTestCase(Helmholtz3DIntEqTestCase): check_gradient = True +class SphereIntEqTestCase(IntEqTestCase): + resolutions = [2, 1] + name = "sphere" + + def get_mesh(self, resolution, target_order): + from meshmode.mesh.generation import generate_icosphere + from meshmode.mesh.refinement import refine_uniformly + mesh = refine_uniformly( + generate_icosphere(1, target_order), + resolution) + + return mesh + + 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 + + class MergedCubesIntEqTestCase(Helmholtz3DIntEqTestCase): resolutions = [1.4] name = "merged-cubes" diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 83f78d64d2291bd3dfa676e3e0c2a77c16409dc0..6894e15e1a5be78324fcae96edabc7f4d32f8420 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()'