diff --git a/examples/exterior-stokes-2d.py b/examples/exterior-stokes-2d.py index 331fb295101ffd8ad9523798899b920d0fc7e9e6..44545bf80a421557b28b49d5b80e0d7cbea8dcdc 100644 --- a/examples/exterior-stokes-2d.py +++ b/examples/exterior-stokes-2d.py @@ -22,6 +22,7 @@ ovsmp_target_order = 4*target_order qbx_order = 4 fmm_order = 7 mu = 1 +circle_rad = 1.5 # }}} @@ -36,7 +37,7 @@ def main(): from meshmode.mesh.generation import ( # noqa make_curve_mesh, starfish, ellipse, drop) mesh = make_curve_mesh( - lambda t: ellipse(1, t), + lambda t: circle_rad * ellipse(1, t), np.linspace(0, 1, nelements+1), target_order) coarse_density_discr = Discretization( @@ -63,9 +64,8 @@ def main(): 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) + meanless_sigma_sym = cse(sigma_sym - sym.mean(2, 1, sigma_sym)) + int_sigma = sym.Ones() * sym.integral(2, 1, sigma_sym) nvec_sym = sym.make_sym_vector("normal", dim) mu_sym = sym.var("mu") @@ -76,10 +76,9 @@ def main(): 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 ) + bdry_op_sym = -loc_sign * 0.5 * sigma_sym - stresslet_obj.apply(sigma_sym, nvec_sym, mu_sym, + qbx_forced_limit='avg') + stokeslet_obj.apply(meanless_sigma_sym, mu_sym, + qbx_forced_limit='avg') - (0.5/np.pi) * int_sigma # }}} @@ -105,8 +104,8 @@ def main(): #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 + xcomp = (-cl.clmath.log(r) + (x - loc[0])**2/r**2) * scaling - (y - loc[1])*strength*0.125/r**2 + 3.3 + ycomp = ((x - loc[0])*(y - loc[1])/r**2) * scaling + (x - loc[0])*strength*0.125/r**2 + 1.5 return [ xcomp, ycomp ] @@ -119,7 +118,7 @@ def main(): 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) + bvp_rhs = bind(qbx, (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), @@ -133,13 +132,13 @@ def main(): # {{{ 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) + int_val = bind(qbx, sym.integral(2, 1, sigma_sym))(queue, sigma=sigma) + int_val = -int_val/(2 * 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 + representation_sym = - stresslet_obj.apply(sigma_sym, nvec_sym, mu_sym, qbx_forced_limit=2) + stokeslet_obj.apply( + meanless_sigma_sym, mu_sym, qbx_forced_limit=2) - u_A_sym_vol + sigma_int_val_sym from sumpy.visualization import FieldPlotter @@ -149,16 +148,16 @@ def main(): 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 circle_mask(test_points, radius): + return (test_points[0,:]**2 + test_points[1,:]**2 > radius**2) - def outside_circle(test_points): - mask = circle_mask(test_points) + def outside_circle(test_points, radius): + mask = circle_mask(test_points, radius) return np.array([ row[mask] for row in test_points]) - eval_points = outside_circle(eval_points) + eval_points = outside_circle(eval_points, radius=circle_rad) from pytential.target import PointsTarget vel = bind( (qbx, PointsTarget(eval_points)), @@ -167,7 +166,7 @@ def main(): from sumpy.visualization import FieldPlotter fplot = FieldPlotter(np.zeros(2), extent=6, npoints=100) - plot_pts = outside_circle(fplot.points) + plot_pts = outside_circle(fplot.points, radius=circle_rad) plot_vel = bind( (qbx, PointsTarget(plot_pts)), representation_sym)(queue, sigma=sigma, mu=mu, normal=normal, sigma_int_val=int_val, omega=omega) @@ -197,7 +196,7 @@ def main(): 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) + mask = circle_mask(fplot.points, radius=circle_rad) for i, vel in enumerate(plot_vel): full_pot[i,mask] = vel.get() diff --git a/pytential/symbolic/stokes.py b/pytential/symbolic/stokes.py index 96c17b501d6ff95a63b1fd7c5347aae041520ed9..adfc23d5eaa5e50144b242e23ada9d6b8e23024d 100644 --- a/pytential/symbolic/stokes.py +++ b/pytential/symbolic/stokes.py @@ -39,7 +39,7 @@ class StokesletWrapper(object): 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, + the vector velocity resulting from convolution with the vector density, and is meant to work similarly to calling S() (which is IntG()). @@ -265,7 +265,7 @@ class StressletWrapper(object): 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 + the kernel with a vector density, and is meant to work similarly to calling S() (which is IntG()). Similar functions are available for other useful things related to