From 4d5675ff9f67f791050e9a4d1da854ad71681692 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Fri, 17 Feb 2017 18:04:34 -0600 Subject: [PATCH 01/13] give option of limiting iterations of refiner when creating qbx object --- pytential/qbx/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index ad60bb22..1fcfc8f1 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,7 @@ 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 -- GitLab From 2110143e5a61fef09ac9ba480b99d85f2938d785 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Fri, 3 Mar 2017 14:51:00 -0600 Subject: [PATCH 02/13] add 2D interior Stokes example --- examples/stokes-2d-interior.py | 220 +++++++++++++++++++++++++++++++++ pytential/symbolic/stokes.py | 174 ++++++++++++++++++++++++++ 2 files changed, 394 insertions(+) create mode 100644 examples/stokes-2d-interior.py create mode 100644 pytential/symbolic/stokes.py diff --git a/examples/stokes-2d-interior.py b/examples/stokes-2d-interior.py new file mode 100644 index 00000000..45ef62cb --- /dev/null +++ b/examples/stokes-2d-interior.py @@ -0,0 +1,220 @@ +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 = 2 + +# 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 = dp * ((y+(h/2.))**2 - h * (y+(h/2.))) + ycomp = 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]]) + + + 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/symbolic/stokes.py b/pytential/symbolic/stokes.py new file mode 100644 index 00000000..00a06f9e --- /dev/null +++ b/pytential/symbolic/stokes.py @@ -0,0 +1,174 @@ +from __future__ import division, absolute_import +import numpy as np + +from pytential import sym +from sumpy.kernel import StokesletKernel, StressletKernel + + +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-indiced 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()). + + .. 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): + """ 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 < 0.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 + + +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-indiced 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()). + + .. 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): + """ Returns an object array of symbolic expressions for the vector resulting from integrating + the dyadic Stresslet kernel with direction vector variable name *dir_vec_sym* and + variable *density_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 < 0.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 + + -- GitLab From 2503530f320e945a99b9a48b9ce4ef107bb8ecc8 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Fri, 3 Mar 2017 15:12:42 -0600 Subject: [PATCH 03/13] flake8 fixes and add exterior stokes example --- examples/exterior-stokes-2d.py | 211 +++++++++++++++++++++++++++++++++ pytential/symbolic/stokes.py | 159 ++++++++++++++----------- 2 files changed, 302 insertions(+), 68 deletions(-) create mode 100644 examples/exterior-stokes-2d.py diff --git a/examples/exterior-stokes-2d.py b/examples/exterior-stokes-2d.py new file mode 100644 index 00000000..01cd03f8 --- /dev/null +++ b/examples/exterior-stokes-2d.py @@ -0,0 +1,211 @@ +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 + +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/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index 00a06f9e..e7939042 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -6,24 +6,25 @@ from sumpy.kernel import StokesletKernel, StressletKernel class StokesletWrapper(object): - """ Wrapper class for the Stokeslet kernel. This class is meant to + """ 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-indiced Stokeslet kernel applied to + in the expansion of the double-indiced Stokeslet kernel applied to the density vector. The object is created - to do some of the set-up and bookkeeping once, rather than every + 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. + 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 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()). + and is meant to work similarly to + calling S() (which is IntG()). - .. 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. + .. 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. """ @@ -31,39 +32,46 @@ class StokesletWrapper(object): 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) } + 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) } + 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") + raise ValueError("unsupported dimension given to StokesletWrapper") def apply(self, density_vec_sym, mu_sym, qbx_forced_limit): - """ Returns an object array of symbolic expressions for the vector resulting from integrating - the dyadic Stokeslet kernel with variable *density_vec_sym*. + """ 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. + :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 + # 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 @@ -73,35 +81,38 @@ class StokesletWrapper(object): ctr_key = tuple(var_ctr) if i < 0.1: - sym_expr[comp] = sym.IntG(self.kernel_dict[ctr_key], density_vec_sym[i], - qbx_forced_limit=qbx_forced_limit, mu=mu_sym) + 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) + 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 class StressletWrapper(object): - """ Wrapper class for the Stresslet kernel. This class is meant to + """ 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-indiced 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 + 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. + 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()). + calling S() (which is IntG()). - .. 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. + .. 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. """ @@ -109,22 +120,26 @@ class StressletWrapper(object): 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) } + 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) } + 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: @@ -132,16 +147,18 @@ class StressletWrapper(object): def apply(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): - """ Returns an object array of symbolic expressions for the vector resulting from integrating - the dyadic Stresslet kernel with direction vector variable name *dir_vec_sym* and + """ Returns an object array of symbolic expressions for the vector + resulting from integrating the dyadic Stresslet kernel with + direction vector variable name *dir_vec_sym* and variable *density_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 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. + :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 @@ -150,7 +167,8 @@ class StressletWrapper(object): for comp in range(self.dim): - # Start variable count for kernel with 1 for the requested result component + # 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 @@ -161,13 +179,18 @@ class StressletWrapper(object): ctr_key = tuple(var_ctr) if i + j < 0.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) + 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) - + 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 -- GitLab From 583ba610d6c5d71532550fe61bd42f1bb4a2085a Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Fri, 3 Mar 2017 15:28:43 -0600 Subject: [PATCH 04/13] flake8 fixes --- pytential/symbolic/stokes.py | 81 ++++++++++++++++-------------------- 1 file changed, 36 insertions(+), 45 deletions(-) diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index e7939042..184d2a47 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -23,8 +23,8 @@ class StokesletWrapper(object): .. 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. + StokesletKernel(icomp=1, jcomp=0) -- and avoid creating multiple expansions + for the same kernel in a different ordering. """ @@ -33,31 +33,29 @@ class StokesletWrapper(object): 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) + (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) + 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): """ 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 @@ -73,7 +71,7 @@ class StokesletWrapper(object): # 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 + base_count[comp] += 1 for i in range(self.dim): var_ctr = base_count.copy() @@ -89,11 +87,9 @@ class StokesletWrapper(object): 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 - 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 @@ -121,30 +117,28 @@ class StressletWrapper(object): 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) + (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) + (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") - + raise ValueError("unsupported dimension given to StressletWrapper") def apply(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): """ Returns an object array of symbolic expressions for the vector @@ -162,7 +156,7 @@ class StressletWrapper(object): """ import itertools - + sym_expr = np.empty((self.dim,), dtype=object) for comp in range(self.dim): @@ -170,7 +164,7 @@ class StressletWrapper(object): # 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 + base_count[comp] += 1 for i, j in itertools.product(range(self.dim), range(self.dim)): var_ctr = base_count.copy() @@ -186,12 +180,9 @@ class StressletWrapper(object): 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) - + 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 - - -- GitLab From 095252485e6674d5b80480522d6c37b08168c54e Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Fri, 3 Mar 2017 15:31:31 -0600 Subject: [PATCH 05/13] add note about exterior representation to comments --- examples/exterior-stokes-2d.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/examples/exterior-stokes-2d.py b/examples/exterior-stokes-2d.py index 01cd03f8..331fb295 100644 --- a/examples/exterior-stokes-2d.py +++ b/examples/exterior-stokes-2d.py @@ -9,6 +9,11 @@ from meshmode.discretization.poly_element import \ 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 -- GitLab From 12a0d2761a91e824a572d55c2516cd57b130b9e3 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Fri, 3 Mar 2017 15:34:20 -0600 Subject: [PATCH 06/13] yet more flake8 fixes --- pytential/symbolic/stokes.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index 184d2a47..990d5365 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -53,9 +53,9 @@ class StokesletWrapper(object): def apply(self, density_vec_sym, mu_sym, qbx_forced_limit): """ Returns an object array of symbolic expressions for the vector - resulting from integrating the dyadic Stokeslet kernel with + 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 @@ -87,9 +87,10 @@ class StokesletWrapper(object): 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 + 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 -- GitLab From aa1110b604dc591c13aca9d414d9d258544b2ef0 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Tue, 7 Mar 2017 00:59:55 -0600 Subject: [PATCH 07/13] add pressure and derivative options to stresslet wrapper class --- pytential/symbolic/stokes.py | 84 ++++++++++++++++++++++++++++++++++-- 1 file changed, 81 insertions(+), 3 deletions(-) diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index 990d5365..07ff4b7c 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -2,7 +2,7 @@ from __future__ import division, absolute_import import numpy as np from pytential import sym -from sumpy.kernel import StokesletKernel, StressletKernel +from sumpy.kernel import StokesletKernel, StressletKernel, LaplaceKernel class StokesletWrapper(object): @@ -78,7 +78,7 @@ class StokesletWrapper(object): var_ctr[i] += 1 ctr_key = tuple(var_ctr) - if i < 0.1: + 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) @@ -173,7 +173,7 @@ class StressletWrapper(object): var_ctr[j] += 1 ctr_key = tuple(var_ctr) - if i + j < 0.1: + if i + j < 1: sym_expr[comp] = sym.IntG( self.kernel_dict[ctr_key], dir_vec_sym[i] * density_vec_sym[j], @@ -187,3 +187,81 @@ class StressletWrapper(object): mu=mu_sym) return sym_expr + + def apply_pressure(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): + """ Returns a symbolic expression to compute the pressure field asssociated + with this stresslet flow. + """ + + 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): + """ Returns an object array of symbolic expressions for the vector + resulting from integrating the derivative of the + dyadic Stresslet kernel with direction vector + variable name *dir_vec_sym* and + 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 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 + 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 -- GitLab From 0293f9e2a9cc1ceda523bba2f2ad6896ca417353 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Tue, 7 Mar 2017 10:43:20 -0600 Subject: [PATCH 08/13] flake8 fixes; add copyright line to stokes.py --- pytential/qbx/__init__.py | 3 ++- pytential/symbolic/stokes.py | 32 ++++++++++++++++++-------------- 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 1fcfc8f1..71a72ce6 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -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), maxiter=maxiter) + InterpolatoryQuadratureSimplexGroupFactory(target_order), + maxiter=maxiter) return lpot, connection @property diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index 07ff4b7c..4d7923af 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -1,4 +1,6 @@ from __future__ import division, absolute_import +__copyright__ = "Copyright (C) 2017 Natalie Beams" + import numpy as np from pytential import sym @@ -8,7 +10,7 @@ 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-indiced Stokeslet kernel applied to + 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, @@ -94,7 +96,7 @@ class StokesletWrapper(object): 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-indiced Stresslet kernel applied to both + 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, @@ -196,29 +198,31 @@ class StressletWrapper(object): 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))) + 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)))) - + 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): + def apply_derivative(self, deriv_dir, density_vec_sym, dir_vec_sym, + mu_sym, qbx_forced_limit): """ Returns an object array of symbolic expressions for the vector resulting from integrating the derivative of the - dyadic Stresslet kernel with direction vector + dyadic Stresslet kernel (wrt target, not source) with direction vector variable name *dir_vec_sym* and variable *density_vec_sym*. -- GitLab From 1f923b362419ce00f0ca32aac1cee4f728b024b1 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Tue, 7 Mar 2017 14:26:50 -0600 Subject: [PATCH 09/13] add pressure computation for Stokeslet --- pytential/symbolic/stokes.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index 4d7923af..69859b36 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -92,6 +92,27 @@ class StokesletWrapper(object): return sym_expr + def apply_pressure(self, density_vec_sym, mu_sym, qbx_forced_limit): + """ Returns a symbolic expression to compute the pressure field asssociated + with this stokeslet flow. + """ + + import itertools + from pytential.symbolic.mappers import DerivativeTaker + kernel = LaplaceKernel(dim=self.dim) + + for i, j in itertools.product(range(self.dim), range(self.dim)): + + if i + j < 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 class StressletWrapper(object): """ Wrapper class for the Stresslet kernel. This class is meant to -- GitLab From 7eaf22e614251bb1402bf366729e95c44478ccc0 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Wed, 8 Mar 2017 00:16:25 -0600 Subject: [PATCH 10/13] add applied stress calculations to stokes wrappers --- pytential/symbolic/stokes.py | 194 ++++++++++++++++++++++++++++++++++- 1 file changed, 189 insertions(+), 5 deletions(-) diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index 69859b36..376dbdb0 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -1,6 +1,26 @@ 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 @@ -22,7 +42,12 @@ class StokesletWrapper(object): and is meant to work similarly to calling S() (which is IntG()). - .. attribute:: kernel_dict + 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 @@ -94,16 +119,16 @@ class StokesletWrapper(object): def apply_pressure(self, density_vec_sym, mu_sym, qbx_forced_limit): """ Returns a symbolic expression to compute the pressure field asssociated - with this stokeslet flow. + with this Stokeslet flow. """ import itertools from pytential.symbolic.mappers import DerivativeTaker kernel = LaplaceKernel(dim=self.dim) - for i, j in itertools.product(range(self.dim), range(self.dim)): + for i in range(self.dim): - if i + j < 1: + if i < 1: sym_expr = DerivativeTaker(i).map_int_g( sym.S(kernel, density_vec_sym[i], qbx_forced_limit=qbx_forced_limit)) @@ -114,6 +139,117 @@ class StokesletWrapper(object): return sym_expr + def apply_derivative(self, deriv_dir, density_vec_sym, + mu_sym, qbx_forced_limit): + """ 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. + """ + + 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 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): + """ 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 @@ -128,6 +264,11 @@ class StressletWrapper(object): 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 @@ -249,7 +390,7 @@ class StressletWrapper(object): :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 direction 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, @@ -290,3 +431,46 @@ class StressletWrapper(object): mu=mu_sym)) return sym_expr + + def apply_stress(self, density_vec_sym, normal_vec_sym, dir_vec_sym, + mu_sym, qbx_forced_limit): + """ 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 -- GitLab From 06f16d4e403793719dce4203523ff2619713720f Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Wed, 8 Mar 2017 00:22:20 -0600 Subject: [PATCH 11/13] flake8 fixes --- pytential/symbolic/stokes.py | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index 376dbdb0..2ffe7131 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -42,9 +42,9 @@ class StokesletWrapper(object): and is meant to work similarly to calling S() (which is IntG()). - Similar functions are available for other useful things related to + 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 + apply_stress (applies symmetric viscous stress tensor in the requested direction). .. attribute:: kernel_dict @@ -122,7 +122,6 @@ class StokesletWrapper(object): with this Stokeslet flow. """ - import itertools from pytential.symbolic.mappers import DerivativeTaker kernel = LaplaceKernel(dim=self.dim) @@ -155,9 +154,8 @@ class StokesletWrapper(object): 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): @@ -192,7 +190,7 @@ class StokesletWrapper(object): def apply_stress(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): """ Returns a vector of symbolic expressions for the force resulting - from the viscous stress: + from the viscous stress: -pressure * I + mu * ( grad U + (grad U).T)), applied in the direction of *dir_vec_sym*. @@ -206,7 +204,7 @@ class StokesletWrapper(object): 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 @@ -220,7 +218,7 @@ class StokesletWrapper(object): 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 @@ -264,9 +262,9 @@ class StressletWrapper(object): 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 + 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 + apply_stress (applies symmetric viscous stress tensor in the requested direction). .. attribute:: kernel_dict @@ -435,7 +433,7 @@ class StressletWrapper(object): def apply_stress(self, density_vec_sym, normal_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): """ Returns a vector of symbolic expressions for the force resulting - from the viscous stress: + from the viscous stress: -pressure * I + mu * ( grad U + (grad U).T)), applied in the direction of *dir_vec_sym*. @@ -455,7 +453,7 @@ class StressletWrapper(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, + 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): @@ -470,7 +468,7 @@ class StressletWrapper(object): sym_expr[comp] = sym_expr[comp] + ( dir_vec_sym[j] * mu_sym * ( sym_grad_matrix[comp][j] + - sym_grad_matrix[j][comp] ) + sym_grad_matrix[j][comp]) ) return sym_expr -- GitLab From a5b9fe2f9f6e2328470f60d58332fb1431f8107e Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Wed, 8 Mar 2017 00:40:25 -0600 Subject: [PATCH 12/13] update interior stokes test to calculate pressure and applied stress --- examples/stokes-2d-interior.py | 32 +++++++++++++++++++++++++++++--- 1 file changed, 29 insertions(+), 3 deletions(-) diff --git a/examples/stokes-2d-interior.py b/examples/stokes-2d-interior.py index 45ef62cb..4f6e7bee 100644 --- a/examples/stokes-2d-interior.py +++ b/examples/stokes-2d-interior.py @@ -20,7 +20,7 @@ target_order = 8 ovsmp_target_order = 4*target_order qbx_order = 2 fmm_order = 7 -mu = 2 +mu = 3 # Test solution type -- either 'fundamental' or 'couette' (default is couette) soln_type = 'couette' @@ -109,8 +109,8 @@ def main(nelements): def couette_soln(x, y, dp, h): scaling = 1./(2*mu) - xcomp = dp * ((y+(h/2.))**2 - h * (y+(h/2.))) - ycomp = 0*y + xcomp = scaling * dp * ((y+(h/2.))**2 - h * (y+(h/2.))) + ycomp = scaling * 0*y return [xcomp, ycomp] @@ -192,6 +192,32 @@ def main(nelements): 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) -- GitLab From 7d7d2dcd236eacabe4562f4ffa85086bcef7d3f6 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Wed, 8 Mar 2017 10:12:58 -0600 Subject: [PATCH 13/13] fix formatting --- pytential/symbolic/stokes.py | 253 ++++++++++++++++++----------------- 1 file changed, 131 insertions(+), 122 deletions(-) diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index 2ffe7131..96c17b50 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -28,31 +28,32 @@ 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. + """ 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()). + 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). + 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. + 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): @@ -79,16 +80,18 @@ class StokesletWrapper(object): raise ValueError("unsupported dimension given to StokesletWrapper") def apply(self, density_vec_sym, mu_sym, qbx_forced_limit): - """ 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. + """ 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) @@ -118,9 +121,7 @@ class StokesletWrapper(object): return sym_expr def apply_pressure(self, density_vec_sym, mu_sym, qbx_forced_limit): - """ Returns a symbolic expression to compute the pressure field asssociated - with this Stokeslet flow. - """ + """ Symbolic expression for pressure field associated with the Stokeslet""" from pytential.symbolic.mappers import DerivativeTaker kernel = LaplaceKernel(dim=self.dim) @@ -140,18 +141,20 @@ class StokesletWrapper(object): def apply_derivative(self, deriv_dir, density_vec_sym, mu_sym, qbx_forced_limit): - """ 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. + """ 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 @@ -189,29 +192,31 @@ class StokesletWrapper(object): def apply_stress(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): - """ 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. + """ 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 @@ -249,25 +254,27 @@ class StokesletWrapper(object): 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. + """ 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()). + 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). + 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 @@ -304,18 +311,19 @@ class StressletWrapper(object): raise ValueError("unsupported dimension given to StressletWrapper") def apply(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): - """ Returns an object array of symbolic expressions for the vector - resulting from integrating the dyadic Stresslet kernel with - direction vector variable name *dir_vec_sym* and - variable *density_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. + """ 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 @@ -351,9 +359,7 @@ class StressletWrapper(object): return sym_expr def apply_pressure(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): - """ Returns a symbolic expression to compute the pressure field asssociated - with this stresslet flow. - """ + """ Symbolic expression for pressure field associated with the stresslet""" import itertools from pytential.symbolic.mappers import DerivativeTaker @@ -380,20 +386,21 @@ class StressletWrapper(object): def apply_derivative(self, deriv_dir, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): - """ Returns an object array of symbolic expressions for the vector - resulting from integrating the derivative of the - dyadic Stresslet kernel (wrt target, not source) with direction vector - variable name *dir_vec_sym* and - 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 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. + """ 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 @@ -432,20 +439,22 @@ class StressletWrapper(object): def apply_stress(self, density_vec_sym, normal_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit): - """ 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. + """ 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) -- GitLab