diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 2ab458a2b12b6d749791b8b23d771a12a627b069..1911a48d1df66c1ea12c7fbb1e8e930ff3a18028 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -523,15 +523,27 @@ def refine_for_global_qbx(lpot_source, wrangler, # {{{ first stage refinement - def visualize_refinement(niter, stage, flags): + def visualize_refinement(niter, stage_nr, stage_name, flags): if not visualize: return - discr = lpot_source.density_discr + if stage_nr == 1: + discr = lpot_source.density_discr + elif stage_nr == 2: + discr = lpot_source.stage2_density_discr + else: + raise ValueError("unexpected stage number") + + flags = flags.get() + logger.info("for stage %s: splitting %d/%d stage-%d elements", + stage_name, np.sum(flags), discr.mesh.nelements, stage_nr) + from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(wrangler.queue, discr, 3) - flags = flags.get().astype(np.bool) + assert len(flags) == discr.mesh.nelements + + flags = flags.astype(np.bool) nodes_flags = np.zeros(discr.nnodes) for grp in discr.groups: meg = grp.mesh_el_group @@ -549,7 +561,7 @@ def refine_for_global_qbx(lpot_source, wrangler, wrangler.queue).as_vector(dtype=object) vis_data.append(("bdry_normals", bdry_normals),) - vis.write_vtk_file("refinement-%03d-%s.vtu" % (niter, stage), vis_data) + vis.write_vtk_file("refinement-%s-%03d.vtu" % (stage_name, niter), vis_data) def warn_max_iterations(): from warnings import warn @@ -602,7 +614,8 @@ def refine_for_global_qbx(lpot_source, wrangler, if violates_kernel_length_scale: iter_violated_criteria.append("kernel length scale") - visualize_refinement(niter, "kernel-length-scale", refine_flags) + visualize_refinement( + niter, 1, "kernel-length-scale", refine_flags) if scaled_max_curvature_threshold is not None: with ProcessLogger(logger, @@ -625,7 +638,7 @@ def refine_for_global_qbx(lpot_source, wrangler, if violates_scaled_max_curv: iter_violated_criteria.append("curvature") - visualize_refinement(niter, "curvature", refine_flags) + visualize_refinement(niter, 1, "curvature", refine_flags) if not iter_violated_criteria: # Only start building trees once the simple length-based criteria @@ -643,7 +656,7 @@ def refine_for_global_qbx(lpot_source, wrangler, refine_flags, debug) if has_disturbed_expansions: iter_violated_criteria.append("disturbed expansions") - visualize_refinement(niter, "disturbed-expansions", refine_flags) + visualize_refinement(niter, 1, "disturbed-expansions", refine_flags) del tree del peer_lists @@ -689,7 +702,7 @@ def refine_for_global_qbx(lpot_source, wrangler, lpot_source, tree, peer_lists, refine_flags, debug) if has_insufficient_quad_res: iter_violated_criteria.append("insufficient quadrature resolution") - visualize_refinement(niter, "quad-resolution", refine_flags) + visualize_refinement(niter, 2, "quad-resolution", refine_flags) if iter_violated_criteria: violated_criteria.append(" and ".join(iter_violated_criteria)) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index ecbdc89926887740e7ccd29554fa6eec491c7bc9..17a1701346b3e199ca80135a8fc490d27d429131 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -648,6 +648,24 @@ def _small_mat_eigenvalues(mat): "eigenvalue formula for %dx%d matrices" % (m, n)) +def _equilateral_parametrization_derivative_matrix(ambient_dim, dim=None, + where=None): + if dim is None: + dim = ambient_dim - 1 + + pder_mat = parametrization_derivative_matrix(ambient_dim, dim, where) + + # The above procedure works well only when the 'reference' end of the + # mapping is in equilateral coordinates. + from modepy.tools import EQUILATERAL_TO_UNIT_MAP + equi_to_unit = EQUILATERAL_TO_UNIT_MAP[dim].a + + # This is the Jacobian of the (equilateral reference element) -> (global) map. + return cse( + np.dot(pder_mat, equi_to_unit), + "equilateral_pder_mat") + + def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, where=None, with_elementwise_max=True): """Return the largest factor by which the reference-to-global @@ -669,17 +687,8 @@ def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, where=None, # stretching, where 'stretching' (*2 to remove bi-unit reference element) # reflects available quadrature resolution in that direction. - pder_mat = parametrization_derivative_matrix(ambient_dim, dim, where) - - # The above procedure works well only when the 'reference' end of the - # mapping is in equilateral coordinates. - from modepy.tools import EQUILATERAL_TO_UNIT_MAP - equi_to_unit = EQUILATERAL_TO_UNIT_MAP[dim].a - - # This is the Jacobian of the (equilateral reference element) -> (global) map. - equi_pder_mat = cse( - np.dot(pder_mat, equi_to_unit), - "equilateral_pder_mat") + equi_pder_mat = _equilateral_parametrization_derivative_matrix( + ambient_dim, dim, where) # Compute eigenvalues of J^T to compute SVD. equi_pder_mat_jtj = cse( diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py index 0d9fca2871c95f7070876e560f9cf94bafaccbc5..a734c074bf115598c54513ba8e15f0a14147561f 100644 --- a/test/test_scalar_int_eq.py +++ b/test/test_scalar_int_eq.py @@ -36,9 +36,11 @@ from functools import partial from meshmode.mesh.generation import ( # noqa ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle, make_curve_mesh) +from meshmode.discretization.visualization import make_visualizer from sumpy.visualization import FieldPlotter from sumpy.symbolic import USE_SYMENGINE from pytential import bind, sym +from pytential.qbx import QBXTargetAssociationFailedException import logging logger = logging.getLogger(__name__) @@ -284,7 +286,106 @@ class ElliptiplaneIntEqTestCase(IntEqTestCase): return perform_flips(mesh, np.ones(mesh.nelements)) inner_radius = 0.2 - outer_radius = 12 + outer_radius = 12 # was '-13' in some large-scale run (?) + + +class BetterplaneIntEqTestCase(IntEqTestCase): + name = "betterplane" + + resolutions = [0.3] + + fmm_backend = "fmmlib" + use_refinement = True + + qbx_order = 3 + fmm_tol = 1e-4 + target_order = 6 + check_gradient = False + check_tangential_deriv = False + + visualize_geometry = True + + #scaled_max_curvature_threshold = 1 + expansion_disturbance_tolerance = 0.3 + + # We're only expecting three digits based on FMM settings. Who are we + # kidding? + gmres_tol = 1e-5 + + def get_mesh(self, resolution, target_order): + from pytools import download_from_web_if_not_present + + download_from_web_if_not_present( + "https://raw.githubusercontent.com/inducer/geometries/a869fc3/" + "surface-3d/betterplane.brep") + + from meshmode.mesh.io import generate_gmsh, ScriptWithFilesSource + mesh = generate_gmsh( + ScriptWithFilesSource(""" + Merge "betterplane.brep"; + + Mesh.CharacteristicLengthMax = %(lcmax)f; + Mesh.ElementOrder = 2; + Mesh.CharacteristicLengthExtendFromBoundary = 0; + + // 2D mesh optimization + // Mesh.Lloyd = 1; + + l_superfine() = Unique(Abs(Boundary{ Surface{ + 27, 25, 17, 13, 18 }; })); + l_fine() = Unique(Abs(Boundary{ Surface{ 2, 6, 7}; })); + l_coarse() = Unique(Abs(Boundary{ Surface{ 14, 16 }; })); + + // p() = Unique(Abs(Boundary{ Line{l_fine()}; })); + // Characteristic Length{p()} = 0.05; + + Field[1] = Attractor; + Field[1].NNodesByEdge = 100; + Field[1].EdgesList = {l_superfine()}; + + Field[2] = Threshold; + Field[2].IField = 1; + Field[2].LcMin = 0.075; + Field[2].LcMax = %(lcmax)f; + Field[2].DistMin = 0.1; + Field[2].DistMax = 0.4; + + Field[3] = Attractor; + Field[3].NNodesByEdge = 100; + Field[3].EdgesList = {l_fine()}; + + Field[4] = Threshold; + Field[4].IField = 3; + Field[4].LcMin = 0.1; + Field[4].LcMax = %(lcmax)f; + Field[4].DistMin = 0.15; + Field[4].DistMax = 0.4; + + Field[5] = Attractor; + Field[5].NNodesByEdge = 100; + Field[5].EdgesList = {l_coarse()}; + + Field[6] = Threshold; + Field[6].IField = 5; + Field[6].LcMin = 0.15; + Field[6].LcMax = %(lcmax)f; + Field[6].DistMin = 0.2; + Field[6].DistMax = 0.4; + + Field[7] = Min; + Field[7].FieldsList = {2, 4, 6}; + + Background Field = 7; + """ % { + "lcmax": resolution, + }, ["betterplane.brep"]), 2) + + # Flip elements--gmsh generates inside-out geometry. + from meshmode.mesh.processing import perform_flips + return perform_flips(mesh, np.ones(mesh.nelements)) + + inner_radius = 0.2 + outer_radius = 15 # }}} @@ -322,19 +423,39 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): qbx_lpot_kwargs["fmm_order"] = case.qbx_order + 5 qbx = QBXLayerPotentialSource( - pre_density_discr, fine_order=source_order, qbx_order=case.qbx_order, + pre_density_discr, + fine_order=source_order, + qbx_order=case.qbx_order, + _from_sep_smaller_min_nsources_cumul=30, fmm_backend=case.fmm_backend, **qbx_lpot_kwargs) if case.use_refinement: if case.k != 0: refiner_extra_kwargs["kernel_length_scale"] = 5/case.k + if hasattr(case, "scaled_max_curvature_threshold"): + refiner_extra_kwargs["_scaled_max_curvature_threshold"] = \ + case.scaled_max_curvature_threshold + + if hasattr(case, "expansion_disturbance_tolerance"): + refiner_extra_kwargs["_expansion_disturbance_tolerance"] = \ + case.expansion_disturbance_tolerance + + if hasattr(case, "refinement_maxiter"): + refiner_extra_kwargs["maxiter"] = case.refinement_maxiter + + #refiner_extra_kwargs["visualize"] = True + print("%d elements before refinement" % pre_density_discr.mesh.nelements) qbx, _ = qbx.with_refinement(**refiner_extra_kwargs) print("%d elements after refinement" % qbx.density_discr.mesh.nelements) density_discr = qbx.density_discr + if hasattr(case, "visualize_geometry") and case.visualize_geometry: + bdry_vis = make_visualizer(queue, density_discr, case.target_order) + bdry_vis.write_vtk_file("geometry.vtu", []) + # {{{ plot geometry if 0: @@ -349,7 +470,6 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): pt.show() elif mesh.ambient_dim == 3: - from meshmode.discretization.visualization import make_visualizer bdry_vis = make_visualizer(queue, density_discr, case.target_order+3) bdry_normals = bind(density_discr, sym.normal(3))(queue)\ @@ -456,13 +576,22 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): rhs = bind(density_discr, op.prepare_rhs(sym.var("bc")))(queue, bc=bc) - from pytential.solve import gmres - gmres_result = gmres( - bound_op.scipy_op(queue, "u", dtype, **concrete_knl_kwargs), - rhs, - tol=case.gmres_tol, - progress=True, - hard_failure=True) + try: + from pytential.solve import gmres + gmres_result = gmres( + bound_op.scipy_op(queue, "u", dtype, **concrete_knl_kwargs), + rhs, + tol=case.gmres_tol, + progress=True, + hard_failure=True, + stall_iterations=50, no_progress_factor=1.05) + except QBXTargetAssociationFailedException as e: + bdry_vis = make_visualizer(queue, density_discr, case.target_order+3) + + bdry_vis.write_vtk_file("failed-targets-%s.vtu" % resolution, [ + ("failed_targets", e.failed_target_flags), + ]) + raise print("gmres state:", gmres_result.state) u = gmres_result.solution @@ -584,10 +713,9 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): # }}} - # {{{ 3D plotting + # {{{ any-D file plotting - if visualize and qbx.ambient_dim == 3: - from meshmode.discretization.visualization import make_visualizer + if visualize: bdry_vis = make_visualizer(queue, density_discr, case.target_order+3) bdry_normals = bind(density_discr, sym.normal(qbx.ambient_dim))(queue)\ @@ -607,7 +735,6 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): qbx_tgt_tol = qbx.copy(target_association_tolerance=0.15) from pytential.target import PointsTarget - from pytential.qbx import QBXTargetAssociationFailedException try: solved_pot = bind( @@ -633,7 +760,6 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): [ ("solved_pot", solved_pot), ("true_pot", true_pot), - ("pot_diff", solved_pot-true_pot), ] )