From f5118f3ea341ac44939d72749b7ea19c6543d937 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 3 Jun 2021 23:03:40 -0500 Subject: [PATCH] Get nodes/normals directly from the discretization collection --- examples/advection/surface.py | 4 +- examples/advection/var-velocity.py | 2 +- examples/advection/weak.py | 4 +- examples/geometry.py | 6 +-- examples/maxwell/cavities.py | 2 +- examples/wave/var-propagation-speed.py | 4 +- examples/wave/wave-min-mpi.py | 2 +- examples/wave/wave-min.py | 2 +- examples/wave/wave-op-mpi.py | 4 +- examples/wave/wave-op-var-velocity.py | 4 +- examples/wave/wave-op.py | 4 +- grudge/models/advection.py | 6 +-- grudge/models/em.py | 4 +- grudge/models/wave.py | 8 +-- test/test_grudge.py | 75 ++++++++------------------ test/test_mpi_communication.py | 4 +- test/test_op.py | 8 +-- 17 files changed, 57 insertions(+), 86 deletions(-) diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 9dcdc93e..c81e2338 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -172,7 +172,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False): # {{{ Surface advection operator # velocity field - x = thaw(op.nodes(dcoll), actx) + x = thaw(dcoll.nodes(), actx) c = make_obj_array([-x[1], x[0], 0.0])[:dim] def f_initial_condition(x): @@ -234,7 +234,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False): df = dof_desc.DOFDesc(FACE_RESTR_INTERIOR) face_discr = dcoll.discr_from_dd(df) - face_normal = thaw(op.normal(dcoll, dd=df), actx) + face_normal = thaw(dcoll.normal(dd=df), actx) from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(actx, face_discr) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index b050c60b..fca97d27 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -169,7 +169,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False): from grudge.models.advection import VariableCoefficientAdvectionOperator - x = thaw(op.nodes(dcoll), actx) + x = thaw(dcoll.nodes(), actx) # velocity field if dim == 1: diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 9300ab01..0d0abb6a 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -152,13 +152,13 @@ def main(ctx_factory, dim=2, order=4, visualize=False): dcoll, c, inflow_u=lambda t: u_analytic( - thaw(op.nodes(dcoll, dd=BTAG_ALL), actx), + thaw(dcoll.nodes(dd=BTAG_ALL), actx), t=t ), flux_type=flux_type ) - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) u = u_analytic(nodes, t=0) def rhs(t, u): diff --git a/examples/geometry.py b/examples/geometry.py index faf0fc42..0c7fb113 100644 --- a/examples/geometry.py +++ b/examples/geometry.py @@ -51,9 +51,9 @@ def main(write_output=True): dcoll = DiscretizationCollection(actx, mesh, order=4) - nodes = thaw(op.nodes(dcoll), actx) - bdry_nodes = thaw(op.nodes(dcoll, dd=BTAG_ALL), actx) - bdry_normals = thaw(op.normal(dcoll, dd=BTAG_ALL), actx) + nodes = thaw(dcoll.nodes(), actx) + bdry_nodes = thaw(dcoll.nodes(dd=BTAG_ALL), actx) + bdry_normals = thaw(dcoll.normal(dd=BTAG_ALL), actx) if write_output: vis = shortcuts.make_visualizer(dcoll) diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index 3d68e513..5791a3d8 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -84,7 +84,7 @@ def main(ctx_factory, dim=3, order=4, visualize=False): else: return get_rectangular_cavity_mode(actx, x, t, 1, (2, 3)) - fields = cavity_mode(thaw(op.nodes(dcoll), actx), t=0) + fields = cavity_mode(thaw(dcoll.nodes(), actx), t=0) maxwell_operator.check_bc_coverage(mesh) diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index 0b8d51b9..04605b7e 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -63,7 +63,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False): source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim] source_width = 0.05 source_omega = 3 - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) source_center_dist = flat_obj_array( [nodes[i] - source_center[i] for i in range(dcoll.dim)] ) @@ -75,7 +75,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False): ) ) - x = thaw(op.nodes(dcoll), actx) + x = thaw(dcoll.nodes(), actx) ones = dcoll.zeros(actx) + 1 c = actx.np.where(np.dot(x, x) < 0.15, 0.1 * ones, 0.2 * ones) diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py index a3e75c7b..c803a897 100644 --- a/examples/wave/wave-min-mpi.py +++ b/examples/wave/wave-min-mpi.py @@ -84,7 +84,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False): source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim] source_width = 0.05 source_omega = 3 - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) source_center_dist = flat_obj_array( [nodes[i] - source_center[i] for i in range(dcoll.dim)] ) diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py index 7b18d3df..ab983681 100644 --- a/examples/wave/wave-min.py +++ b/examples/wave/wave-min.py @@ -65,7 +65,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False): source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim] source_width = 0.05 source_omega = 3 - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) source_center_dist = flat_obj_array( [nodes[i] - source_center[i] for i in range(dcoll.dim)] ) diff --git a/examples/wave/wave-op-mpi.py b/examples/wave/wave-op-mpi.py index 3fb49d7d..cc8578c9 100644 --- a/examples/wave/wave-op-mpi.py +++ b/examples/wave/wave-op-mpi.py @@ -54,7 +54,7 @@ def wave_flux(dcoll, c, w_tpair): u = w_tpair[0] v = w_tpair[1:] - normal = thaw(op.normal(dcoll, w_tpair.dd), u.int.array_context) + normal = thaw(dcoll.normal(w_tpair.dd), u.int.array_context) flux_weak = flat_obj_array( np.dot(v.avg, normal), @@ -136,7 +136,7 @@ def bump(actx, dcoll, t=0): source_width = 0.05 source_omega = 3 - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) center_dist = flat_obj_array([ nodes[i] - source_center[i] for i in range(dcoll.dim) diff --git a/examples/wave/wave-op-var-velocity.py b/examples/wave/wave-op-var-velocity.py index e48c8d71..07b40b8d 100644 --- a/examples/wave/wave-op-var-velocity.py +++ b/examples/wave/wave-op-var-velocity.py @@ -56,7 +56,7 @@ def wave_flux(dcoll, c, w_tpair): u = w_tpair[0] v = w_tpair[1:] - normal = thaw(op.normal(dcoll, dd), u.int.array_context) + normal = thaw(dcoll.normal(dd), u.int.array_context) flux_weak = flat_obj_array( np.dot(v.avg, normal), @@ -146,7 +146,7 @@ def bump(actx, dcoll, t=0, width=0.05, center=None): center = center[:dcoll.dim] source_omega = 3 - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) center_dist = flat_obj_array([ nodes[i] - center[i] for i in range(dcoll.dim) diff --git a/examples/wave/wave-op.py b/examples/wave/wave-op.py index c0494d9f..56a90171 100644 --- a/examples/wave/wave-op.py +++ b/examples/wave/wave-op.py @@ -52,7 +52,7 @@ def wave_flux(dcoll, c, w_tpair): u = w_tpair[0] v = w_tpair[1:] - normal = thaw(op.normal(dcoll, w_tpair.dd), u.int.array_context) + normal = thaw(dcoll.normal(w_tpair.dd), u.int.array_context) flux_weak = flat_obj_array( np.dot(v.avg, normal), @@ -127,7 +127,7 @@ def bump(actx, dcoll, t=0): source_width = 0.05 source_omega = 3 - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) center_dist = flat_obj_array([ nodes[i] - source_center[i] for i in range(dcoll.dim) diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 2469498a..0ba4acb3 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -43,7 +43,7 @@ def advection_weak_flux(dcoll, flux_type, u_tpair, velocity): """ actx = u_tpair.int.array_context dd = u_tpair.dd - normal = thaw(op.normal(dcoll, dd), actx) + normal = thaw(dcoll.normal(dd), actx) v_dot_n = np.dot(velocity, normal) flux_type = flux_type.lower() @@ -92,7 +92,7 @@ class StrongAdvectionOperator(AdvectionOperatorBase): def flux(self, u_tpair): actx = u_tpair.int.array_context dd = u_tpair.dd - normal = thaw(op.normal(self.dcoll, dd), actx) + normal = thaw(self.dcoll.normal(dd), actx) v_dot_normal = np.dot(self.v, normal) return u_tpair.int * v_dot_normal - self.weak_flux(u_tpair) @@ -285,7 +285,7 @@ def v_dot_n_tpair(actx, dcoll, velocity, trace_dd): from grudge.trace_pair import TracePair from meshmode.discretization.connection import FACE_RESTR_INTERIOR - normal = thaw(op.normal(dcoll, trace_dd.with_discr_tag(None)), actx) + normal = thaw(dcoll.normal(trace_dd.with_discr_tag(None)), actx) v_dot_n = velocity.dot(normal) i = op.project(dcoll, trace_dd.with_discr_tag(None), trace_dd, v_dot_n) diff --git a/grudge/models/em.py b/grudge/models/em.py index 38b73d3d..8960e0d3 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -120,7 +120,7 @@ class MaxwellOperator(HyperbolicOperator): As per Hesthaven and Warburton page 433. """ - normal = thaw(op.normal(self.dcoll, wtpair.dd), self.dcoll._setup_actx) + normal = thaw(self.dcoll.normal(wtpair.dd), self.dcoll._setup_actx) if self.fixed_material: e, h = self.split_eh(wtpair) @@ -220,7 +220,7 @@ class MaxwellOperator(HyperbolicOperator): absorbing boundary conditions. """ - absorb_normal = thaw(op.normal(self.dcoll, dd=self.absorb_tag), + absorb_normal = thaw(self.dcoll.normal(dd=self.absorb_tag), self.dcoll._setup_actx) e, h = self.split_eh(w) diff --git a/grudge/models/wave.py b/grudge/models/wave.py index e101fd04..fa428796 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -90,7 +90,7 @@ class WeakWaveOperator(HyperbolicOperator): u = wtpair[0] v = wtpair[1:] actx = u.int.array_context - normal = thaw(op.normal(self.dcoll, wtpair.dd), actx) + normal = thaw(self.dcoll.normal(wtpair.dd), actx) central_flux_weak = -self.c*flat_obj_array( np.dot(v.avg, normal), @@ -133,7 +133,7 @@ class WeakWaveOperator(HyperbolicOperator): neu_bc = flat_obj_array(neu_u, -neu_v) # radiation BCs ------------------------------------------------------- - rad_normal = thaw(op.normal(dcoll, dd=self.radiation_tag), actx) + rad_normal = thaw(dcoll.normal(dd=self.radiation_tag), actx) rad_u = op.project(dcoll, "vol", self.radiation_tag, u) rad_v = op.project(dcoll, "vol", self.radiation_tag, v) @@ -233,7 +233,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): u = wtpair[1] v = wtpair[2:] actx = u.int.array_context - normal = thaw(op.normal(self.dcoll, wtpair.dd), actx) + normal = thaw(self.dcoll.normal(wtpair.dd), actx) flux_central_weak = -0.5 * flat_obj_array( np.dot(v.int*c.int + v.ext*c.ext, normal), @@ -282,7 +282,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator): neu_bc = flat_obj_array(neu_c, neu_u, -neu_v) # radiation BCs ------------------------------------------------------- - rad_normal = thaw(op.normal(dcoll, dd=self.radiation_tag), actx) + rad_normal = thaw(dcoll.normal(dd=self.radiation_tag), actx) rad_c = op.project(dcoll, "vol", self.radiation_tag, self.c) rad_u = op.project(dcoll, "vol", self.radiation_tag, u) diff --git a/test/test_grudge.py b/test/test_grudge.py index 4963642e..522028fe 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -250,7 +250,9 @@ def test_mass_surface_area(actx_factory, name): # }}} # compute max element size - h_max = op.h_max_from_volume(dcoll) + from grudge.dt_utils import h_max_from_volume + + h_max = h_max_from_volume(dcoll) eoc.add_data_point(h_max, area_error) @@ -321,7 +323,9 @@ def test_mass_operator_inverse(actx_factory, name): # }}} # compute max element size - h_max = op.h_max_from_volume(dcoll) + from grudge.dt_utils import h_max_from_volume + + h_max = h_max_from_volume(dcoll) eoc.add_data_point(h_max, inv_error) @@ -378,7 +382,7 @@ def test_face_normal_surface(actx_factory, mesh_name): ) surf_normal = surf_normal / actx.np.sqrt(sum(surf_normal**2)) - face_normal_i = thaw(op.normal(dcoll, df), actx) + face_normal_i = thaw(dcoll.normal(df), actx) face_normal_e = dcoll.opposite_face_connection()(face_normal_i) if mesh.ambient_dim == 3: @@ -498,7 +502,7 @@ def test_2d_gauss_theorem(actx_factory): int_1 = op.integral(dcoll, "vol", op.local_div(dcoll, f_volm)) prj_f = op.project(dcoll, "vol", BTAG_ALL, f_volm) - normal = thaw(op.normal(dcoll, BTAG_ALL), actx) + normal = thaw(dcoll.normal(BTAG_ALL), actx) int_2 = op.integral(dcoll, BTAG_ALL, prj_f.dot(normal)) assert abs(int_1 - int_2) < 1e-13 @@ -600,14 +604,14 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False): ambient_dim = dcoll.ambient_dim # variables - f_num = f(thaw(op.nodes(dcoll, dd=dd), actx)) - f_quad_num = f(thaw(op.nodes(dcoll, dd=dq), actx)) + f_num = f(thaw(dcoll.nodes(dd=dd), actx)) + f_quad_num = f(thaw(dcoll.nodes(dd=dq), actx)) from grudge.geometry import normal, summed_curvature kappa = summed_curvature(actx, dcoll, dd=dq) normal = normal(actx, dcoll, dd=dq) - face_normal = thaw(op.normal(dcoll, df), actx) + face_normal = thaw(dcoll.normal(df), actx) face_f = op.project(dcoll, dd, df, f_num) # operators @@ -627,7 +631,9 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False): logger.info("errors: global %.5e local %.5e", err_global, err_local) # compute max element size - h_max = op.h_max_from_volume(dcoll) + from grudge.dt_utils import h_max_from_volume + + h_max = h_max_from_volume(dcoll) eoc_global.add_data_point(h_max, err_global) eoc_local.add_data_point(h_max, err_local) @@ -747,12 +753,12 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ "weak": WeakAdvectionOperator}[op_type] adv_operator = op_class(dcoll, v, inflow_u=lambda t: u_analytic( - thaw(op.nodes(dcoll, dd=BTAG_ALL), actx), + thaw(dcoll.nodes(dd=BTAG_ALL), actx), t=t ), flux_type=flux_type) - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) u = u_analytic(nodes, t=0) def rhs(t, u): @@ -763,7 +769,9 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ else: final_time = 0.2 - h_max = op.h_max_from_volume(dcoll, dim=dcoll.ambient_dim) + from grudge.dt_utils import h_max_from_volume + + h_max = h_max_from_volume(dcoll, dim=dcoll.ambient_dim) dt = dt_factor * h_max/order**2 nsteps = (final_time // dt) + 1 dt = final_time/nsteps + 1e-15 @@ -842,7 +850,7 @@ def test_convergence_maxwell(actx_factory, order): def analytic_sol(x, t=0): return get_rectangular_cavity_mode(actx, x, t, 1, (1, 2, 2)) - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) fields = analytic_sol(nodes, t=0) from grudge.models.em import MaxwellOperator @@ -939,7 +947,7 @@ def test_improvement_quadrature(actx_factory, order): discr_tag_to_group_factory=discr_tag_to_group_factory ) - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) def zero_inflow(dtag, t=0): dd = dof_desc.DOFDesc(dtag, qtag) @@ -988,7 +996,7 @@ def test_bessel(actx_factory): dcoll = DiscretizationCollection(actx, mesh, order=3) - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) r = actx.np.sqrt(nodes[0]**2 + nodes[1]**2) # FIXME: Bessel functions need to brought out of the symbolic @@ -1079,50 +1087,13 @@ def test_empty_boundary(actx_factory): a=(-0.5,)*dim, b=(0.5,)*dim, nelements_per_axis=(8,)*dim, order=4) dcoll = DiscretizationCollection(actx, mesh, order=4) - normal = op.normal(dcoll, BTAG_NONE) + normal = dcoll.normal(BTAG_NONE) from meshmode.dof_array import DOFArray for component in normal: assert isinstance(component, DOFArray) assert len(component) == len(dcoll.discr_from_dd(BTAG_NONE).groups) -# {{{ DiscretizationCollection testing - -def test_dcoll_nodes_and_normals(actx_factory): - actx = actx_factory() - - from meshmode.mesh import BTAG_ALL - from meshmode.mesh.generation import generate_warped_rect_mesh - - mesh = generate_warped_rect_mesh(dim=2, order=4, nelements_side=6) - - dcoll = DiscretizationCollection(actx, mesh, order=3) - - # Nodes should be *indentical* - nodes = thaw(op.nodes(dcoll), actx) - dcoll_nodes = thaw(dcoll.nodes(), actx) - - assert op.norm(dcoll, nodes - dcoll_nodes, np.inf) == 0.0 - - nodes_btag = thaw(op.nodes(dcoll, BTAG_ALL), actx) - dcoll_nodes_btag = thaw(dcoll.nodes(BTAG_ALL), actx) - - assert op.norm(dcoll, nodes_btag - dcoll_nodes_btag, np.inf) == 0.0 - - # Normals should be *indentical* - normals = thaw(op.normal(dcoll, "int_faces"), actx) - dcoll_normals = thaw(dcoll.normal("int_faces"), actx) - - assert op.norm(dcoll, normals - dcoll_normals, np.inf) == 0.0 - - normals_btag = thaw(op.normal(dcoll, BTAG_ALL), actx) - dcoll_normals_btag = thaw(dcoll.normal(BTAG_ALL), actx) - - assert op.norm(dcoll, normals_btag - dcoll_normals_btag, np.inf) == 0.0 - -# }}} - - # You can test individual routines by typing # $ python test_grudge.py 'test_routine()' diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py index 879d3bf8..4e3484e9 100644 --- a/test/test_mpi_communication.py +++ b/test/test_mpi_communication.py @@ -76,7 +76,7 @@ def simple_mpi_communication_entrypoint(): dcoll = DiscretizationCollection(actx, local_mesh, order=5, mpi_communicator=comm) - x = thaw(op.nodes(dcoll), actx) + x = thaw(dcoll.nodes(), actx) myfunc = actx.np.sin(np.dot(x, [2, 3])) from grudge.dof_desc import as_dofdesc @@ -147,7 +147,7 @@ def mpi_communication_entrypoint(): source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim] source_width = 0.05 source_omega = 3 - nodes = thaw(op.nodes(dcoll), actx) + nodes = thaw(dcoll.nodes(), actx) source_center_dist = flat_obj_array( [nodes[i] - source_center[i] for i in range(dcoll.dim)] ) diff --git a/test/test_op.py b/test/test_op.py index 2d12409e..6aa3ba79 100644 --- a/test/test_op.py +++ b/test/test_op.py @@ -88,7 +88,7 @@ def test_gradient(actx_factory, form, dim, order, vectorize, nested, result[dim-1] = result[dim-1] * (-np.pi/2*actx.np.sin(np.pi/2*x[dim-1])) return result - x = thaw(op.nodes(dcoll), actx) + x = thaw(dcoll.nodes(), actx) if vectorize: u = make_obj_array([(i+1)*f(x) for i in range(dim)]) @@ -98,7 +98,7 @@ def test_gradient(actx_factory, form, dim, order, vectorize, nested, def get_flux(u_tpair): dd = u_tpair.dd dd_allfaces = dd.with_dtag("all_faces") - normal = thaw(op.normal(dcoll, dd), actx) + normal = thaw(dcoll.normal(dd), actx) u_avg = u_tpair.avg if vectorize: if nested: @@ -212,7 +212,7 @@ def test_divergence(actx_factory, form, dim, order, vectorize, nested, result = result + deriv return result - x = thaw(op.nodes(dcoll), actx) + x = thaw(dcoll.nodes(), actx) if vectorize: u = make_obj_array([(i+1)*f(x) for i in range(dim)]) @@ -224,7 +224,7 @@ def test_divergence(actx_factory, form, dim, order, vectorize, nested, def get_flux(u_tpair): dd = u_tpair.dd dd_allfaces = dd.with_dtag("all_faces") - normal = thaw(op.normal(dcoll, dd), actx) + normal = thaw(dcoll.normal(dd), actx) flux = u_tpair.avg @ normal return op.project(dcoll, dd, dd_allfaces, flux) -- GitLab