diff --git a/examples/exterior-stokes-2d.py b/examples/exterior-stokes-2d.py new file mode 100644 index 0000000000000000000000000000000000000000..331fb295101ffd8ad9523798899b920d0fc7e9e6 --- /dev/null +++ b/examples/exterior-stokes-2d.py @@ -0,0 +1,216 @@ +import numpy as np +import pyopencl as cl +import pyopencl.clmath # noqa + +from meshmode.discretization import Discretization +from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + +from pytential import bind, sym, norm # noqa + +# {{{ set some constants for use below +# This program tests an exterior Stokes flow in 2D using the +# compound representation given in Hsiao & Kress, +# ``On an integral equation for the two-dimensional exterior Stokes problem,'' +# Applied Numerical Mathematics 1 (1985). + + +nelements = 50 +mesh_order = 4 +target_order = 8 +ovsmp_target_order = 4*target_order +qbx_order = 4 +fmm_order = 7 +mu = 1 + +# }}} + + +def main(): + import logging + logging.basicConfig(level=logging.INFO) + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + from meshmode.mesh.generation import ( # noqa + make_curve_mesh, starfish, ellipse, drop) + mesh = make_curve_mesh( + lambda t: ellipse(1, t), + np.linspace(0, 1, nelements+1), + target_order) + coarse_density_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + from pytential.qbx import QBXLayerPotentialSource + stick_out = 0.05 + qbx, _ = QBXLayerPotentialSource( + coarse_density_discr, fine_order=ovsmp_target_order, qbx_order=qbx_order, + fmm_order=fmm_order, target_stick_out_factor=stick_out + ).with_refinement() + + density_discr = qbx.density_discr + normal = bind(density_discr, sym.normal(2).as_vector())(queue) + path_length = bind(density_discr, sym.integral(2, 1, 1))(queue) + + # {{{ describe bvp + + from pytential.symbolic.stokes import StressletWrapper, StokesletWrapper + from pytools.obj_array import make_obj_array + dim=2 + cse = sym.cse + + sigma_sym = sym.make_sym_vector("sigma", dim) + sqrt_w = sym.sqrt_jac_q_weight(2) + inv_sqrt_w_sigma = cse(sigma_sym/sqrt_w) + meanless_inv_sqrt_w_sigma = cse(inv_sqrt_w_sigma - sym.mean(2, 1, inv_sqrt_w_sigma)) + int_inv_sqrt_w_sigma = sym.Ones() * sym.integral(2, 1, inv_sqrt_w_sigma) + + nvec_sym = sym.make_sym_vector("normal", dim) + mu_sym = sym.var("mu") + + # -1 for interior Dirichlet + # +1 for exterior Dirichlet + loc_sign = 1 + + stresslet_obj = StressletWrapper(dim=2) + stokeslet_obj = StokesletWrapper(dim=2) + bdry_op_sym = -loc_sign * 0.5 * sigma_sym - sqrt_w * ( + stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit='avg') + + stokeslet_obj.apply(meanless_inv_sqrt_w_sigma, mu_sym, qbx_forced_limit='avg') + + (0.5/np.pi) * int_inv_sqrt_w_sigma ) + + # }}} + + bound_op = bind(qbx, bdry_op_sym) + + # {{{ fix rhs and solve + + def fund_soln(x, y, loc, strength): + #with direction (1,0) for point source + r = cl.clmath.sqrt((x - loc[0])**2 + (y - loc[1])**2) + scaling = strength/(4*np.pi*mu) + xcomp = (-cl.clmath.log(r) + (x - loc[0])**2/r**2) * scaling + ycomp = ((x - loc[0])*(y - loc[1])/r**2) * scaling + return [ xcomp, ycomp ] + + def rotlet_soln(x, y, loc): + r = cl.clmath.sqrt((x - loc[0])**2 + (y - loc[1])**2) + xcomp = -(y - loc[1])/r**2 + ycomp = (x - loc[0])/r**2 + return [ xcomp, ycomp ] + + def fund_and_rot_soln(x, y, loc, strength): + #with direction (1,0) for point source + r = cl.clmath.sqrt((x - loc[0])**2 + (y - loc[1])**2) + scaling = strength/(4*np.pi*mu) + xcomp = (-cl.clmath.log(r) + (x - loc[0])**2/r**2) * scaling - (y - loc[1])*strength*0.125/r**2 + ycomp = ((x - loc[0])*(y - loc[1])/r**2) * scaling + (x - loc[0])*strength*0.125/r**2 + return [ xcomp, ycomp ] + + + nodes = density_discr.nodes().with_queue(queue) + fund_soln_loc = np.array([0.5, -0.2]) + strength = 100. + bc = fund_and_rot_soln(nodes[0], nodes[1], fund_soln_loc, strength) + + omega_sym = sym.make_sym_vector("omega", dim) + u_A_sym_bdry = stokeslet_obj.apply(omega_sym, mu_sym, qbx_forced_limit=1) + + omega = [cl.array.to_device(queue, (strength/path_length)*np.ones(len(nodes[0]))), cl.array.to_device(queue, np.zeros(len(nodes[0])))] + bvp_rhs = bind(qbx, sqrt_w*(sym.make_sym_vector("bc",dim) - u_A_sym_bdry))(queue, bc=bc, mu=mu, omega=omega) + from pytential.solve import gmres + gmres_result = gmres( + bound_op.scipy_op(queue, "sigma", np.float64, mu=mu, normal=normal), + bvp_rhs, tol=1e-9, progress=True, + stall_iterations=0, + hard_failure=True) + + # }}} + + + # {{{ postprocess/visualize + sigma = gmres_result.solution + sigma_int_val_sym = sym.make_sym_vector("sigma_int_val", 2) + int_val = bind(qbx, sym.integral(2, 1, inv_sqrt_w_sigma))(queue, sigma=sigma) + int_val = -int_val/(np.pi) + print("int_val = ", int_val) + + u_A_sym_vol = stokeslet_obj.apply(omega_sym, mu_sym, qbx_forced_limit=2) + representation_sym = - stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit=2) - stokeslet_obj.apply( + meanless_inv_sqrt_w_sigma, mu_sym, qbx_forced_limit=2) + sigma_int_val_sym + u_A_sym_vol + + + from sumpy.visualization import FieldPlotter + nsamp = 30 + eval_points_1d = np.linspace(-3., 3., nsamp) + eval_points = np.zeros((2, len(eval_points_1d)**2)) + eval_points[0,:] = np.tile(eval_points_1d, len(eval_points_1d)) + eval_points[1,:] = np.repeat(eval_points_1d, len(eval_points_1d)) + + def circle_mask(test_points): + return (test_points[0,:]**2 + test_points[1,:]**2 > 1) + + def outside_circle(test_points): + mask = circle_mask(test_points) + return np.array([ + row[mask] + for row in test_points]) + + eval_points = outside_circle(eval_points) + from pytential.target import PointsTarget + vel = bind( + (qbx, PointsTarget(eval_points)), + representation_sym)(queue, sigma=sigma, mu=mu, normal=normal, sigma_int_val=int_val, omega=omega) + print("@@@@@@@@") + + from sumpy.visualization import FieldPlotter + fplot = FieldPlotter(np.zeros(2), extent=6, npoints=100) + plot_pts = outside_circle(fplot.points) + plot_vel = bind( + (qbx, PointsTarget(plot_pts)), + representation_sym)(queue, sigma=sigma, mu=mu, normal=normal, sigma_int_val=int_val, omega=omega) + + def get_obj_array(obj_array): + from pytools.obj_array import make_obj_array + return make_obj_array([ + ary.get() + for ary in obj_array + ]) + exact_soln = fund_and_rot_soln(cl.array.to_device(queue, eval_points[0]), cl.array.to_device( + queue, eval_points[1]), fund_soln_loc, strength) + vel = get_obj_array(vel) + err = vel-get_obj_array(exact_soln) + rel_err = err/(get_obj_array(exact_soln)) + + print("@@@@@@@@") + print("vel[0], err[0], rel_err[0] ***** vel[1], err[1], rel_err[1]: ") + for i in range(len(vel[0])): + print("%15.8e, %15.8e, %15.8e ***** %15.8e, %15.8e, %15.8e\n"%( + vel[0][i], err[0][i], rel_err[0][i], vel[1][i], err[1][i], rel_err[1][i])) + + print("@@@@@@@@") + print("L2 error estimate: ", np.sqrt((6./(nsamp-1))**2*np.sum(err[0]*err[0]) + (6./(nsamp-1))**2*np.sum(err[1]*err[1]))) + print("L2 rel error estimate: ", np.sqrt((6./(nsamp-1))**2*np.sum(rel_err[0]*rel_err[0]) + (6./(nsamp-1))**2*np.sum(rel_err[1]*rel_err[1]))) + print("max error at sampled points: ", max(abs(err[0])), max(abs(err[1]))) + print("max rel error at sampled points: ", max(abs(rel_err[0])), max(abs(rel_err[1]))) + + full_pot = np.zeros_like(fplot.points) * float("nan") + mask = circle_mask(fplot.points) + + for i, vel in enumerate(plot_vel): + full_pot[i,mask] = vel.get() + + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + plt.quiver(fplot.points[0], fplot.points[1], full_pot[0], full_pot[1], linewidth=0.1) + plt.savefig("exterior-2d-field.pdf") + + # }}} + + + +if __name__ == "__main__": + main() diff --git a/examples/stokes-2d-interior.py b/examples/stokes-2d-interior.py new file mode 100644 index 0000000000000000000000000000000000000000..4f6e7bee4158fc39d0a9f00ad9f276c22604e905 --- /dev/null +++ b/examples/stokes-2d-interior.py @@ -0,0 +1,246 @@ +import numpy as np +import pyopencl as cl +import pyopencl.clmath # noqa + +from meshmode.discretization import Discretization +from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + +from pytential.target import PointsTarget +import matplotlib +matplotlib.use("Agg") + + +from pytential import bind, sym, norm # noqa + +# {{{ set some constants for use below + +mesh_order = 4 +target_order = 8 +ovsmp_target_order = 4*target_order +qbx_order = 2 +fmm_order = 7 +mu = 3 + +# Test solution type -- either 'fundamental' or 'couette' (default is couette) +soln_type = 'couette' + +# }}} + + +def main(nelements): + import logging + logging.basicConfig(level=logging.INFO) + + def get_obj_array(obj_array): + from pytools.obj_array import make_obj_array + return make_obj_array([ + ary.get() + for ary in obj_array + ]) + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + from meshmode.mesh.generation import ( # noqa + make_curve_mesh, starfish, ellipse, drop) + mesh = make_curve_mesh( + lambda t: starfish(t), + np.linspace(0, 1, nelements+1), + target_order) + coarse_density_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + from pytential.qbx import QBXLayerPotentialSource + stick_out = 0.05 + qbx, _ = QBXLayerPotentialSource( + coarse_density_discr, fine_order=ovsmp_target_order, qbx_order=qbx_order, + fmm_order=fmm_order, target_stick_out_factor=stick_out + ).with_refinement() + + density_discr = qbx.density_discr + nodes = density_discr.nodes().with_queue(queue) + + # Get normal vectors for the density discretization -- used in integration with stresslet + mv_normal = bind(density_discr, sym.normal(2))(queue) + normal = mv_normal.as_vector(np.object) + + + # {{{ describe bvp + + from sumpy.kernel import LaplaceKernel + from pytential.symbolic.stokes import StressletWrapper + from pytools.obj_array import make_obj_array + dim=2 + cse = sym.cse + + nvec_sym = sym.make_sym_vector("normal", dim) + sigma_sym = sym.make_sym_vector("sigma", dim) + mu_sym = sym.var("mu") + sqrt_w = sym.sqrt_jac_q_weight(2) + inv_sqrt_w_sigma = cse(sigma_sym/sqrt_w) + + # -1 for interior Dirichlet + # +1 for exterior Dirichlet + loc_sign = -1 + + # Create stresslet object + stresslet_obj = StressletWrapper(dim=2) + + # Describe boundary operator + bdry_op_sym = loc_sign * 0.5 * sigma_sym + sqrt_w * stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit='avg') + + # Bind to the qbx discretization + bound_op = bind(qbx, bdry_op_sym) + + # }}} + + + # {{{ fix rhs and solve + + def fund_soln(x, y, loc): + #with direction (1,0) for point source + r = cl.clmath.sqrt((x - loc[0])**2 + (y - loc[1])**2) + scaling = 1./(4*np.pi*mu) + xcomp = (-cl.clmath.log(r) + (x - loc[0])**2/r**2) * scaling + ycomp = ((x - loc[0])*(y - loc[1])/r**2) * scaling + return [ xcomp, ycomp ] + + def couette_soln(x, y, dp, h): + scaling = 1./(2*mu) + xcomp = scaling * dp * ((y+(h/2.))**2 - h * (y+(h/2.))) + ycomp = scaling * 0*y + return [xcomp, ycomp] + + + + if soln_type == 'fundamental': + pt_loc = np.array([2.0, 0.0]) + bc = fund_soln(nodes[0], nodes[1], pt_loc) + else: + dp = -10. + h = 2.5 + bc = couette_soln(nodes[0], nodes[1], dp, h) + + # Get rhs vector + bvp_rhs = bind(qbx, sqrt_w*sym.make_sym_vector("bc",dim))(queue, bc=bc) + + from pytential.solve import gmres + gmres_result = gmres( + bound_op.scipy_op(queue, "sigma", np.float64, mu=mu, normal=normal), + bvp_rhs, tol=1e-9, progress=True, + stall_iterations=0, + hard_failure=True) + + # }}} + + # {{{ postprocess/visualize + sigma = gmres_result.solution + + # Describe representation of solution for evaluation in domain + representation_sym = stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit=-2) + + from sumpy.visualization import FieldPlotter + nsamp = 10 + eval_points_1d = np.linspace(-1., 1., nsamp) + eval_points = np.zeros((2, len(eval_points_1d)**2)) + eval_points[0,:] = np.tile(eval_points_1d, len(eval_points_1d)) + eval_points[1,:] = np.repeat(eval_points_1d, len(eval_points_1d)) + + + gamma_sym = sym.var("gamma") + inv_sqrt_w_gamma = cse(gamma_sym/sqrt_w) + constant_laplace_rep = sym.D(LaplaceKernel(dim=2), inv_sqrt_w_gamma, qbx_forced_limit=None) + sqrt_w_vec = bind(qbx, sqrt_w)(queue) + + def general_mask(test_points): + const_density = bind((qbx, PointsTarget(test_points)), constant_laplace_rep)(queue, gamma=sqrt_w_vec).get() + return (abs(const_density) > 0.1) + + def inside_domain(test_points): + mask = general_mask(test_points) + return np.array([ + row[mask] + for row in test_points]) + + def stride_hack(arr): + from numpy.lib.stride_tricks import as_strided + return np.array(as_strided(arr, strides=(8 * len(arr[0]), 8))) + + eval_points = inside_domain(eval_points) + eval_points_dev = cl.array.to_device(queue, eval_points) + + # Evaluate the solution at the evaluation points + vel = bind( + (qbx, PointsTarget(eval_points_dev)), + representation_sym)(queue, sigma=sigma, mu=mu, normal=normal) + print("@@@@@@@@") + vel = get_obj_array(vel) + + + if soln_type == 'fundamental': + exact_soln = fund_soln(eval_points_dev[0], eval_points_dev[1], pt_loc) + else: + exact_soln = couette_soln(eval_points_dev[0], eval_points_dev[1], dp, h) + err = vel - get_obj_array(exact_soln) + + print("@@@@@@@@") + + print("L2 error estimate: ", np.sqrt((2./(nsamp-1))**2*np.sum(err[0]*err[0]) + (2./(nsamp-1))**2*np.sum(err[1]*err[1]))) + max_error_loc = [abs(err[0]).argmax(), abs(err[1]).argmax()] + print("max error at sampled points: ", max(abs(err[0])), max(abs(err[1]))) + print("exact velocity at max error points: x -> ", err[0][max_error_loc[0]], ", y -> ", err[1][max_error_loc[1]]) + + from pytential.symbolic.mappers import DerivativeTaker + rep_pressure = stresslet_obj.apply_pressure(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit=-2) + pressure = bind((qbx, PointsTarget(eval_points_dev)), + rep_pressure)(queue, sigma=sigma, mu=mu, normal=normal) + pressure = pressure.get() + print "pressure = ", pressure + + x_dir_vecs = np.zeros((2,len(eval_points[0]))) + x_dir_vecs[0,:] = 1.0 + y_dir_vecs = np.zeros((2, len(eval_points[0]))) + y_dir_vecs[1,:] = 1.0 + x_dir_vecs = cl.array.to_device(queue, x_dir_vecs) + y_dir_vecs = cl.array.to_device(queue, y_dir_vecs) + dir_vec_sym = sym.make_sym_vector("force_direction", dim) + rep_stress = stresslet_obj.apply_stress(inv_sqrt_w_sigma, nvec_sym, dir_vec_sym, mu_sym, qbx_forced_limit=-2) + + applied_stress_x = bind((qbx, PointsTarget(eval_points_dev)), + rep_stress)(queue, sigma=sigma, normal=normal, force_direction=x_dir_vecs, mu=mu) + applied_stress_x = get_obj_array(applied_stress_x) + applied_stress_y = bind((qbx, PointsTarget(eval_points_dev)), + rep_stress)(queue, sigma=sigma, normal=normal, force_direction=y_dir_vecs, mu=mu) + applied_stress_y = get_obj_array(applied_stress_y) + + print "stress applied to x direction: ", applied_stress_x + print "stress applied to y direction: ", applied_stress_y + + + import matplotlib.pyplot as plt + plt.quiver(eval_points[0], eval_points[1], vel[0], vel[1], linewidth=0.1) + file_name = "field-n%s.pdf"%(nelements) + plt.savefig(file_name) + + return (max(abs(err[0])), max(abs(err[1]))) + + # }}} + + +if __name__ == "__main__": + n_elements = np.array([30, 60, 120]) + max_errs_x = np.zeros(3) + max_errs_y = np.zeros(3) + + + for i in range(3): + max_errs_x[i], max_errs_y[i] = main(n_elements[i]) + + + print("@@@@@@@@@@@@@@@@@@@@@") + print("max errs in x for each mesh: ", max_errs_x) + print("max errs in y for each mesh: ", max_errs_y) + + diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index ad60bb22c2930fe987102e78e22143e7674817a0..71a72ce6cccfe34101a9ce245234942251627677 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -248,7 +248,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): + (group.nelements,)) @memoize_method - def with_refinement(self, target_order=None): + def with_refinement(self, target_order=None, maxiter=3): """ :returns: a tuple ``(lpot_src, cnx)``, where ``lpot_src`` is a :class:`QBXLayerPotentialSource` and ``cnx`` is a @@ -262,7 +262,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): if target_order is None: target_order = self.density_discr.groups[0].order lpot, connection = refiner(self, - InterpolatoryQuadratureSimplexGroupFactory(target_order)) + InterpolatoryQuadratureSimplexGroupFactory(target_order), + maxiter=maxiter) return lpot, connection @property diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py new file mode 100644 index 0000000000000000000000000000000000000000..96c17b501d6ff95a63b1fd7c5347aae041520ed9 --- /dev/null +++ b/pytential/symbolic/stokes.py @@ -0,0 +1,483 @@ +from __future__ import division, absolute_import +__copyright__ = "Copyright (C) 2017 Natalie Beams" + +__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 + +from pytential import sym +from sumpy.kernel import StokesletKernel, StressletKernel, LaplaceKernel + + +class StokesletWrapper(object): + """ Wrapper class for the Stokeslet kernel. + + This class is meant to shield the user from the messiness of writing + out every term in the expansion of the double-indexed Stokeslet kernel + applied to the density vector. The object is created + to do some of the set-up and bookkeeping once, rather than every + time we want to create a symbolic expression based on the kernel -- say, + once when we solve for the density, and once when we want a symbolic + representation for the solution, for example. + + The apply() function returns the integral expressions needed for + the vector velocity resulting from convolution with the vectory density, + and is meant to work similarly to + calling S() (which is IntG()). + + Similar functions are available for other useful things related to + the flow: apply_pressure, apply_derivative (target derivative), + apply_stress (applies symmetric viscous stress tensor in + the requested direction). + + .. attribute:: kernel_dict + + The dictionary allows us to exploit symmetry -- that + StokesletKernel(icomp=0, jcomp=1) is identical to + StokesletKernel(icomp=1, jcomp=0) -- and avoid creating multiple expansions + for the same kernel in a different ordering. + """ + + def __init__(self, dim=None): + + self.dim = dim + if dim == 2: + self.kernel_dict = { + (2, 0): StokesletKernel(dim=2, icomp=0, jcomp=0), + (1, 1): StokesletKernel(dim=2, icomp=0, jcomp=1), + (0, 2): StokesletKernel(dim=2, icomp=1, jcomp=1) + } + + elif dim == 3: + self.kernel_dict = { + (2, 0, 0): StokesletKernel(dim=3, icomp=0, jcomp=0), + (1, 1, 0): StokesletKernel(dim=3, icomp=0, jcomp=1), + (1, 0, 1): StokesletKernel(dim=3, icomp=0, jcomp=2), + (0, 2, 0): StokesletKernel(dim=3, icomp=1, jcomp=1), + (0, 1, 1): StokesletKernel(dim=3, icomp=1, jcomp=2), + (0, 0, 2): StokesletKernel(dim=3, icomp=2, jcomp=2) + } + + else: + raise ValueError("unsupported dimension given to StokesletWrapper") + + def apply(self, density_vec_sym, mu_sym, qbx_forced_limit): + """ Symbolic expressions for integrating Stokeslet kernel + + Returns an object array of symbolic expressions for the vector + resulting from integrating the dyadic Stokeslet kernel with + variable *density_vec_sym*. + + :arg density_vec_sym: a symbolic vector variable for the density vector + :arg mu_sym: a symbolic variable for the viscosity + :arg qbx_forced_limit: the qbx_forced_limit argument to be passed on + to IntG. +/-1 for exterior/interior one-sided boundary limit, + +/-2 for exterior/interior off-boundary evaluation, and 'avg' + for the average of the two one-sided boundary limits. + """ + + sym_expr = np.empty((self.dim,), dtype=object) + + for comp in range(self.dim): + + # Start variable count for kernel with 1 for the requested result + # component + base_count = np.zeros(self.dim, dtype=np.int) + base_count[comp] += 1 + + for i in range(self.dim): + var_ctr = base_count.copy() + var_ctr[i] += 1 + ctr_key = tuple(var_ctr) + + if i < 1: + sym_expr[comp] = sym.IntG( + self.kernel_dict[ctr_key], density_vec_sym[i], + qbx_forced_limit=qbx_forced_limit, mu=mu_sym) + + else: + sym_expr[comp] = sym_expr[comp] + sym.IntG( + self.kernel_dict[ctr_key], density_vec_sym[i], + qbx_forced_limit=qbx_forced_limit, mu=mu_sym) + + return sym_expr + + def apply_pressure(self, density_vec_sym, mu_sym, qbx_forced_limit): + """ Symbolic expression for pressure field associated with the Stokeslet""" + + from pytential.symbolic.mappers import DerivativeTaker + kernel = LaplaceKernel(dim=self.dim) + + for i in range(self.dim): + + if i < 1: + sym_expr = DerivativeTaker(i).map_int_g( + sym.S(kernel, density_vec_sym[i], + qbx_forced_limit=qbx_forced_limit)) + else: + sym_expr = sym_expr + (DerivativeTaker(i).map_int_g( + sym.S(kernel, density_vec_sym[i], + qbx_forced_limit=qbx_forced_limit))) + + return sym_expr + + def apply_derivative(self, deriv_dir, density_vec_sym, + mu_sym, qbx_forced_limit): + """ Symbolic derivative of velocity from Stokeslet. + + Returns an object array of symbolic expressions for the vector + resulting from integrating the *deriv_dir* derivative of the + dyadic Stokeslet kernel (wrt target, not source) with + variable *density_vec_sym*. + + :arg deriv_dir: which derivative we want: 0, 1, or 2 for x, y, z + :arg density_vec_sym: a symbolic vector variable for the density vector + :arg mu_sym: a symbolic variable for the viscosity + :arg qbx_forced_limit: the qbx_forced_limit argument to be passed + on to IntG. +/-1 for exterior/interior one-sided boundary limit, + +/-2 for exterior/interior off-boundary evaluation, and 'avg' + for the average of the two one-sided boundary limits. + """ + + from pytential.symbolic.mappers import DerivativeTaker + + sym_expr = np.empty((self.dim,), dtype=object) + + for comp in range(self.dim): + + # Start variable count for kernel with 1 for the requested result + # component + base_count = np.zeros(self.dim, dtype=np.int) + base_count[comp] += 1 + + for i in range(self.dim): + var_ctr = base_count.copy() + var_ctr[i] += 1 + ctr_key = tuple(var_ctr) + + if i < 1: + sym_expr[comp] = DerivativeTaker(deriv_dir).map_int_g( + sym.IntG(self.kernel_dict[ctr_key], + density_vec_sym[i], + qbx_forced_limit=qbx_forced_limit, + mu=mu_sym)) + + else: + sym_expr[comp] = sym_expr[comp] + DerivativeTaker( + deriv_dir).map_int_g( + sym.IntG(self.kernel_dict[ctr_key], + density_vec_sym[i], + qbx_forced_limit=qbx_forced_limit, + mu=mu_sym)) + + return sym_expr + + def apply_stress(self, density_vec_sym, dir_vec_sym, + mu_sym, qbx_forced_limit): + """ Symbolic expression for viscous stress applied to direction + + Returns a vector of symbolic expressions for the force resulting + from the viscous stress: + -pressure * I + mu * ( grad U + (grad U).T)), + applied in the direction of *dir_vec_sym*. + + Note that this computation is very similar to computing + a double-layer potential with the stresslet kernel. + The difference is that here the direction vector is the + direction applied to the stress tensor and is applied + outside of the integration, whereas the stresslet calculation + uses the normal vectors at every source point. As such, the + length of the argument passed in for the stresslet velocity + calculation (after binding) is the same length as the number + of source points/nodes; when calling this routine, the number + of direction vectors should be the same as the number of targets. + + :arg density_vec_sym: a symbolic vector variable for the density vector + :arg dir_vec_sym: a symbolic vector for the application direction + :arg mu_sym: a symbolic variable for the viscosity + :arg qbx_forced_limit: the qbx_forced_limit argument to be passed + on to IntG. +/-1 for exterior/interior one-sided boundary limit, + +/-2 for exterior/interior off-boundary evaluation, and 'avg' + for the average of the two one-sided boundary limits. + """ + + import itertools + + sym_expr = np.empty((self.dim,), dtype=object) + stresslet_obj = StressletWrapper(dim=self.dim) + + for comp in range(self.dim): + + # Start variable count for kernel with 1 for the requested result + # component + base_count = np.zeros(self.dim, dtype=np.int) + base_count[comp] += 1 + + for i, j in itertools.product(range(self.dim), range(self.dim)): + var_ctr = base_count.copy() + var_ctr[i] += 1 + var_ctr[j] += 1 + ctr_key = tuple(var_ctr) + + if i + j < 1: + sym_expr[comp] = dir_vec_sym[i] * sym.IntG( + stresslet_obj.kernel_dict[ctr_key], + density_vec_sym[j], + qbx_forced_limit=qbx_forced_limit, mu=mu_sym) + + else: + sym_expr[comp] = sym_expr[comp] + dir_vec_sym[i] * sym.IntG( + stresslet_obj.kernel_dict[ctr_key], + density_vec_sym[j], + qbx_forced_limit=qbx_forced_limit, + mu=mu_sym) + + return sym_expr + + +class StressletWrapper(object): + """ Wrapper class for the Stresslet kernel. + + This class is meant to shield the user from the messiness of writing + out every term in the expansion of the triple-indexed Stresslet + kernel applied to both a normal vector and the density vector. + The object is created to do some of the set-up and bookkeeping once, + rather than every time we want to create a symbolic expression based + on the kernel -- say, once when we solve for the density, and once when + we want a symbolic representation for the solution, for example. + + The apply() function returns the integral expressions needed for convolving + the kernel with a vectory density, and is meant to work similarly to + calling S() (which is IntG()). + + Similar functions are available for other useful things related to + the flow: apply_pressure, apply_derivative (target derivative), + apply_stress (applies symmetric viscous stress tensor in + the requested direction). + + .. attribute:: kernel_dict + + The dictionary allows us to exploit symmetry -- that + StressletKernel(icomp=0, jcomp=1, kcomp=2) is identical to + StressletKernel(icomp=1, jcomp=2, kcomp=0) -- and avoid creating + multiple expansions for the same kernel in a different ordering. + + """ + + def __init__(self, dim=None): + + self.dim = dim + if dim == 2: + self.kernel_dict = { + (3, 0): StressletKernel(dim=2, icomp=0, jcomp=0, kcomp=0), + (2, 1): StressletKernel(dim=2, icomp=0, jcomp=0, kcomp=1), + (1, 2): StressletKernel(dim=2, icomp=0, jcomp=1, kcomp=1), + (0, 3): StressletKernel(dim=2, icomp=1, jcomp=1, kcomp=1) + } + + elif dim == 3: + self.kernel_dict = { + (3, 0, 0): StressletKernel(dim=3, icomp=0, jcomp=0, kcomp=0), + (2, 1, 0): StressletKernel(dim=3, icomp=0, jcomp=0, kcomp=1), + (2, 0, 1): StressletKernel(dim=3, icomp=0, jcomp=0, kcomp=2), + (1, 2, 0): StressletKernel(dim=3, icomp=0, jcomp=1, kcomp=1), + (1, 1, 1): StressletKernel(dim=3, icomp=0, jcomp=1, kcomp=2), + (1, 0, 2): StressletKernel(dim=3, icomp=0, jcomp=2, kcomp=2), + (0, 3, 0): StressletKernel(dim=3, icomp=1, jcomp=1, kcomp=1), + (0, 2, 1): StressletKernel(dim=3, icomp=1, jcomp=1, kcomp=2), + (0, 1, 2): StressletKernel(dim=3, icomp=1, jcomp=2, kcomp=2), + (0, 0, 3): StressletKernel(dim=3, icomp=2, jcomp=2, kcomp=2) + } + + else: + raise ValueError("unsupported dimension given to StressletWrapper") + + def apply(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): + """ Symbolic expressions for integrating stresslet kernel + + Returns an object array of symbolic expressions for the vector + resulting from integrating the dyadic stresslet kernel with + variable *density_vec_sym* and source direction vectors *dir_vec_sym*. + + :arg density_vec_sym: a symbolic vector variable for the density vector + :arg dir_vec_sym: a symbolic vector variable for the direction vector + :arg mu_sym: a symbolic variable for the viscosity + :arg qbx_forced_limit: the qbx_forced_limit argument to be passed + on to IntG. +/-1 for exterior/interior one-sided boundary limit, + +/-2 for exterior/interior off-boundary evaluation, and 'avg' + for the average of the two one-sided boundary limits. + """ + + import itertools + + sym_expr = np.empty((self.dim,), dtype=object) + + for comp in range(self.dim): + + # Start variable count for kernel with 1 for the requested result + # component + base_count = np.zeros(self.dim, dtype=np.int) + base_count[comp] += 1 + + for i, j in itertools.product(range(self.dim), range(self.dim)): + var_ctr = base_count.copy() + var_ctr[i] += 1 + var_ctr[j] += 1 + ctr_key = tuple(var_ctr) + + if i + j < 1: + sym_expr[comp] = sym.IntG( + self.kernel_dict[ctr_key], + dir_vec_sym[i] * density_vec_sym[j], + qbx_forced_limit=qbx_forced_limit, mu=mu_sym) + + else: + sym_expr[comp] = sym_expr[comp] + sym.IntG( + self.kernel_dict[ctr_key], + dir_vec_sym[i] * density_vec_sym[j], + qbx_forced_limit=qbx_forced_limit, + mu=mu_sym) + + return sym_expr + + def apply_pressure(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): + """ Symbolic expression for pressure field associated with the stresslet""" + + import itertools + from pytential.symbolic.mappers import DerivativeTaker + kernel = LaplaceKernel(dim=self.dim) + + factor = (2. * mu_sym) + + for i, j in itertools.product(range(self.dim), range(self.dim)): + + if i + j < 1: + sym_expr = factor * DerivativeTaker(i).map_int_g( + DerivativeTaker(j).map_int_g( + sym.S(kernel, density_vec_sym[i] * dir_vec_sym[j], + qbx_forced_limit=qbx_forced_limit))) + else: + sym_expr = sym_expr + ( + factor * DerivativeTaker(i).map_int_g( + DerivativeTaker(j).map_int_g( + sym.S(kernel, + density_vec_sym[i] * dir_vec_sym[j], + qbx_forced_limit=qbx_forced_limit)))) + + return sym_expr + + def apply_derivative(self, deriv_dir, density_vec_sym, dir_vec_sym, + mu_sym, qbx_forced_limit): + """ Symbolic derivative of velocity from stresslet. + + Returns an object array of symbolic expressions for the vector + resulting from integrating the *deriv_dir* derivative of the + dyadic stresslet kernel (wrt target, not source) with + variable *density_vec_sym* and source direction vectors *dir_vec_sym*. + + :arg deriv_dir: which derivative we want: 0, 1, or 2 for x, y, z + :arg density_vec_sym: a symbolic vector variable for the density vector + :arg dir_vec_sym: a symbolic vector variable for the normal direction + :arg mu_sym: a symbolic variable for the viscosity + :arg qbx_forced_limit: the qbx_forced_limit argument to be passed + on to IntG. +/-1 for exterior/interior one-sided boundary limit, + +/-2 for exterior/interior off-boundary evaluation, and 'avg' + for the average of the two one-sided boundary limits. + """ + + import itertools + from pytential.symbolic.mappers import DerivativeTaker + + sym_expr = np.empty((self.dim,), dtype=object) + + for comp in range(self.dim): + + # Start variable count for kernel with 1 for the requested result + # component + base_count = np.zeros(self.dim, dtype=np.int) + base_count[comp] += 1 + + for i, j in itertools.product(range(self.dim), range(self.dim)): + var_ctr = base_count.copy() + var_ctr[i] += 1 + var_ctr[j] += 1 + ctr_key = tuple(var_ctr) + + if i + j < 1: + sym_expr[comp] = DerivativeTaker(deriv_dir).map_int_g( + sym.IntG(self.kernel_dict[ctr_key], + dir_vec_sym[i] * density_vec_sym[j], + qbx_forced_limit=qbx_forced_limit, mu=mu_sym)) + + else: + sym_expr[comp] = sym_expr[comp] + DerivativeTaker( + deriv_dir).map_int_g( + sym.IntG(self.kernel_dict[ctr_key], + dir_vec_sym[i] * density_vec_sym[j], + qbx_forced_limit=qbx_forced_limit, + mu=mu_sym)) + + return sym_expr + + def apply_stress(self, density_vec_sym, normal_vec_sym, dir_vec_sym, + mu_sym, qbx_forced_limit): + """ Symbolic expression for viscous stress applied to direction + + Returns a vector of symbolic expressions for the force resulting + from the viscous stress: + -pressure * I + mu * ( grad U + (grad U).T)), + applied in the direction of *dir_vec_sym*. + + :arg density_vec_sym: a symbolic vector variable for the density vector + :arg normal_vec_sym: a symbolic vector variable for the normal vectors + (outward facing normals at source locations) + :arg dir_vec_sym: a symbolic vector for the application direction + :arg mu_sym: a symbolic variable for the viscosity + :arg qbx_forced_limit: the qbx_forced_limit argument to be passed + on to IntG. +/-1 for exterior/interior one-sided boundary limit, + +/-2 for exterior/interior off-boundary evaluation, and 'avg' + for the average of the two one-sided boundary limits. + """ + + sym_expr = np.empty((self.dim,), dtype=object) + + # Build velocity derivative matrix + sym_grad_matrix = np.empty((self.dim, self.dim), dtype=object) + for i in range(self.dim): + sym_grad_matrix[:, i] = self.apply_derivative(i, density_vec_sym, + normal_vec_sym, mu_sym, qbx_forced_limit) + + for comp in range(self.dim): + + # First, add the pressure term: + sym_expr[comp] = - dir_vec_sym[comp] * self.apply_pressure( + density_vec_sym, normal_vec_sym, + mu_sym, qbx_forced_limit) + + # Now add the velocity derivative components + for j in range(self.dim): + sym_expr[comp] = sym_expr[comp] + ( + dir_vec_sym[j] * mu_sym * ( + sym_grad_matrix[comp][j] + + sym_grad_matrix[j][comp]) + ) + + return sym_expr