From 6a4411040d0a2ccaf5547b8bdcdf53ec900f2206 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 23 Mar 2017 22:48:51 -0500 Subject: [PATCH 1/3] Add exterior 2D Stokes test --- .../test_stokes.py | 170 ++++++++++++------ 1 file changed, 111 insertions(+), 59 deletions(-) rename examples/exterior-stokes-2d.py => test/test_stokes.py (57%) diff --git a/examples/exterior-stokes-2d.py b/test/test_stokes.py similarity index 57% rename from examples/exterior-stokes-2d.py rename to test/test_stokes.py index 44545bf8..9c602f47 100644 --- a/examples/exterior-stokes-2d.py +++ b/test/test_stokes.py @@ -1,3 +1,28 @@ +from __future__ import division, absolute_import, print_function + +__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 import pyopencl as cl import pyopencl.clmath # noqa @@ -6,34 +31,29 @@ 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). +from pytools.obj_array import make_obj_array +from sumpy.visualization import FieldPlotter +from pytential import bind, sym, norm # noqa +import logging -nelements = 50 -mesh_order = 4 -target_order = 8 -ovsmp_target_order = 4*target_order -qbx_order = 4 -fmm_order = 7 -mu = 1 -circle_rad = 1.5 -# }}} +def run_exterior_stokes_2d(ctx_factory, nelements, + mesh_order=4, target_order=4, qbx_order=4, + fmm_order=10, mu=1, circle_rad=1.5, do_plot=False): + # 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). -def main(): - import logging logging.basicConfig(level=logging.INFO) cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) + ovsmp_target_order = 4*target_order + from meshmode.mesh.generation import ( # noqa make_curve_mesh, starfish, ellipse, drop) mesh = make_curve_mesh( @@ -58,15 +78,13 @@ def main(): # {{{ describe bvp from pytential.symbolic.stokes import StressletWrapper, StokesletWrapper - from pytools.obj_array import make_obj_array - dim=2 + dim = 2 cse = sym.cse sigma_sym = sym.make_sym_vector("sigma", dim) - sqrt_w = sym.sqrt_jac_q_weight(2) 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") @@ -78,7 +96,7 @@ def main(): stokeslet_obj = StokesletWrapper(dim=2) 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 + qbx_forced_limit='avg') - (0.5/np.pi) * int_sigma # }}} @@ -92,13 +110,13 @@ def main(): 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 ] + 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 ] + return [xcomp, ycomp] def fund_and_rot_soln(x, y, loc, strength): #with direction (1,0) for point source @@ -106,14 +124,13 @@ def main(): 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 + 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 ] - + 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) @@ -128,8 +145,8 @@ def main(): # }}} + # {{{ postprocess/visualize - # {{{ 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, sigma_sym))(queue, sigma=sigma) @@ -140,16 +157,14 @@ def main(): 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 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)) + 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, radius): - return (test_points[0,:]**2 + test_points[1,:]**2 > radius**2) + return (test_points[0, :]**2 + test_points[1, :]**2 > radius**2) def outside_circle(test_points, radius): mask = circle_mask(test_points, radius) @@ -161,55 +176,92 @@ def main(): 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) + 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, 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) + 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) + + 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) + + # FIXME: Pointwise relative errors don't make sense! 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])) + if 0: + 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("@@@@@@@@") + + l2_err = np.sqrt((6./(nsamp-1))**2*np.sum(err[0]*err[0]) + (6./(nsamp-1))**2*np.sum(err[1]*err[1])) + l2_rel_err = 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("L2 error estimate: ", l2_err) + print("L2 rel error estimate: ", l2_rel_err) 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, radius=circle_rad) + if do_plot: + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt - for i, vel in enumerate(plot_vel): - full_pot[i,mask] = vel.get() + full_pot = np.zeros_like(fplot.points) * float("nan") + mask = circle_mask(fplot.points, radius=circle_rad) - 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") + for i, vel in enumerate(plot_vel): + full_pot[i, mask] = vel.get() - # }}} + plt.quiver( + fplot.points[0], fplot.points[1], + full_pot[0], full_pot[1], + linewidth=0.1) + plt.savefig("exterior-2d-field.pdf") + + return qbx.h_max, l2_err + + +def test_exterior_stokes_2d(ctx_factory, qbx_order=3): + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + for nelements in [20, 50]: + h_max, l2_err = run_exterior_stokes_2d(ctx_factory, nelements) + eoc_rec.add_data_point(h_max, l2_err) + + print(eoc_rec) + assert eoc_rec.order_estimate() > qbx_order +# You can test individual routines by typing +# $ python test_layer_pot.py 'test_routine()' if __name__ == "__main__": - main() + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from py.test.cmdline import main + main([__file__]) + +# vim: fdm=marker -- GitLab From 597822942d6ff281d54a1a39352a16876565a6e6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 23 Mar 2017 22:53:40 -0500 Subject: [PATCH 2/3] PEP8 Stokes test --- test/test_stokes.py | 54 +++++++++++++++++++++++++++++++-------------- 1 file changed, 38 insertions(+), 16 deletions(-) diff --git a/test/test_stokes.py b/test/test_stokes.py index 9c602f47..b9523fee 100644 --- a/test/test_stokes.py +++ b/test/test_stokes.py @@ -35,6 +35,7 @@ from pytools.obj_array import make_obj_array from sumpy.visualization import FieldPlotter from pytential import bind, sym, norm # noqa +from pytential.solve import gmres import logging @@ -94,9 +95,12 @@ def run_exterior_stokes_2d(ctx_factory, nelements, stresslet_obj = StressletWrapper(dim=2) stokeslet_obj = StokesletWrapper(dim=2) - 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 + 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) # }}} @@ -122,8 +126,12 @@ def run_exterior_stokes_2d(ctx_factory, nelements, #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 + 3.3 - ycomp = ((x - loc[0])*(y - loc[1])/r**2) * scaling + (x - loc[0])*strength*0.125/r**2 + 1.5 + 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] nodes = density_discr.nodes().with_queue(queue) @@ -132,11 +140,14 @@ def run_exterior_stokes_2d(ctx_factory, nelements, 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, (sym.make_sym_vector("bc",dim) + u_A_sym_bdry))(queue, bc=bc, mu=mu, omega=omega) - from pytential.solve import gmres + u_A_sym_bdry = stokeslet_obj.apply(omega_sym, mu_sym, qbx_forced_limit=1) # noqa: N806, E501 + + 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, sym.make_sym_vector("bc", dim) + u_A_sym_bdry + )(queue, bc=bc, mu=mu, omega=omega) gmres_result = gmres( bound_op.scipy_op(queue, "sigma", np.float64, mu=mu, normal=normal), bvp_rhs, tol=1e-9, progress=True, @@ -153,9 +164,13 @@ def run_exterior_stokes_2d(ctx_factory, nelements, 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(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 + u_A_sym_vol = stokeslet_obj.apply(omega_sym, mu_sym, qbx_forced_limit=2) # noqa: N806, E501 + 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) nsamp = 30 eval_points_1d = np.linspace(-3., 3., nsamp) @@ -213,13 +228,18 @@ def run_exterior_stokes_2d(ctx_factory, nelements, print("@@@@@@@@") - l2_err = np.sqrt((6./(nsamp-1))**2*np.sum(err[0]*err[0]) + (6./(nsamp-1))**2*np.sum(err[1]*err[1])) - l2_rel_err = 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])) + l2_err = np.sqrt( + (6./(nsamp-1))**2*np.sum(err[0]*err[0]) + + (6./(nsamp-1))**2*np.sum(err[1]*err[1])) + l2_rel_err = 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("L2 error estimate: ", l2_err) print("L2 rel error estimate: ", l2_rel_err) 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]))) + print("max rel error at sampled points: ", + max(abs(rel_err[0])), max(abs(rel_err[1]))) if do_plot: import matplotlib @@ -238,6 +258,8 @@ def run_exterior_stokes_2d(ctx_factory, nelements, linewidth=0.1) plt.savefig("exterior-2d-field.pdf") + # }}} + return qbx.h_max, l2_err -- GitLab From 4a798dfe43e96be493e0802f61db10512d5ae655 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 23 Mar 2017 23:34:04 -0500 Subject: [PATCH 3/3] Stokes test: Add missing PyOpenCL fixture --- test/test_stokes.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_stokes.py b/test/test_stokes.py index b9523fee..02a71b1b 100644 --- a/test/test_stokes.py +++ b/test/test_stokes.py @@ -34,6 +34,9 @@ from meshmode.discretization.poly_element import \ from pytools.obj_array import make_obj_array from sumpy.visualization import FieldPlotter +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + from pytential import bind, sym, norm # noqa from pytential.solve import gmres import logging -- GitLab