From c14939c3e9d22b8e31044f649ff8d848d0359e7f Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sat, 24 Feb 2018 02:22:00 -0600 Subject: [PATCH 1/6] Implemented WADG for everything, all the time --- grudge/symbolic/mappers/__init__.py | 18 ++++--- test/test_grudge.py | 76 +++++++++++++++++++++++++++++ 2 files changed, 87 insertions(+), 7 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index ca11cbf0..d5b2582e 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -515,19 +515,23 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): if with_jacobian: rec_field = jac_tag * rec_field return sum( - sym.inverse_metric_derivative( - rst_axis, expr.op.xyz_axis, - ambient_dim=self.ambient_dim, dim=self.dim) * - ref_class(rst_axis, dd_in=dd_in)(rec_field) - for rst_axis in range(self.dim)) + ref_class(rst_axis, dd_in=dd_in)( + sym.inverse_metric_derivative( + rst_axis, expr.op.xyz_axis, + ambient_dim=self.ambient_dim, dim=self.dim) *rec_field) + for rst_axis in range(self.dim)) + if isinstance(expr.op, op.MassOperator): return op.RefMassOperator(expr.op.dd_in, expr.op.dd_out)( jac_in * self.rec(expr.field)) elif isinstance(expr.op, op.InverseMassOperator): - return op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_out)( - 1/jac_in * self.rec(expr.field)) + in_inv_op = op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_in) + trans_mass_op = op.RefMassOperator(expr.op.dd_in, expr.op.dd_out) + out_inv_op = op.RefInverseMassOperator(expr.op.dd_out, expr.op.dd_out) + return out_inv_op(trans_mass_op( in_inv_op(self.rec(expr.field))/jac_in )) + elif isinstance(expr.op, op.FaceMassOperator): jac_in_surf = sym.area_element(self.ambient_dim, self.dim - 1, diff --git a/test/test_grudge.py b/test/test_grudge.py index 34532bfc..94456bb7 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -332,6 +332,82 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type assert eoc_rec.order_estimate() > order +@pytest.mark.parametrize("order", [3, 4, 5]) +def test_convergence_curvi(ctx_factory, order): + """Test whether a curvilinear mesh converges""" + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + dims = 2 + ns = [4, 6, 8] + for n in ns: + + from meshmode.mesh.generation import generate_warped_rect_mesh + mesh = generate_warped_rect_mesh(dims, order, n) + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + + advec_v = np.array([0.1, 0.1]) + norm_v = la.norm(advec_v) + + def f(x): + return sym.sin(3*x) + + def u_analytic(x): + return f(-np.dot(advec_v, x)/norm_v+sym.var("t", sym.DD_SCALAR)*norm_v) + + from grudge.models.advection import WeakAdvectionOperator + analytic_sol = bind(discr, u_analytic(sym.nodes(dims))) + + analytic_init = bind(discr, u_analytic(sym.nodes(dims)) ) + + fields = analytic_init(queue, t=0) + + op = WeakAdvectionOperator(advec_v, + inflow_u=u_analytic(sym.nodes(dims, sym.BTAG_ALL) ), + flux_type="central") + + bound_op = bind(discr, op.sym_operator()) + + def rhs(t, u): + return bound_op(queue, t=t, u=u) + + final_t = 0.1 + + h = 1.0/n + dt_factor = 4 + dt = dt_factor * h/order**2 + nsteps = (final_t // dt) + 1 + dt = final_t/nsteps + 1e-15 + + from grudge.shortcuts import set_up_rk4 + dt_stepper = set_up_rk4("u", dt, fields, rhs) + + norm = bind(discr, sym.norm(2, sym.var("u"))) + + step = 0 + for event in dt_stepper.run(t_end=final_t): + if isinstance(event, dt_stepper.StateComputed): + assert event.component_id == "u" + esc = event.state_component + step += 1 + + + sol = analytic_sol(queue, t=step * dt) + total_error = norm(queue, u=(esc - sol)) / norm(queue, u=sol) + eoc_rec.add_data_point(1.0/n, total_error) + + + + print(eoc_rec.pretty_print(abscissa_label="h", + error_label="L2 Error")) + + assert eoc_rec.order_estimate() > order + @pytest.mark.parametrize("order", [3, 4, 5]) def test_convergence_maxwell(ctx_factory, order): -- GitLab From 3ef75d4f1074de825660180aa0ad857f32a6ab38 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sat, 24 Feb 2018 02:28:53 -0600 Subject: [PATCH 2/6] Placate Flake for WADG --- grudge/symbolic/mappers/__init__.py | 10 ++++------ test/test_grudge.py | 8 +++----- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index d5b2582e..e1a26f01 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -516,12 +516,11 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): rec_field = jac_tag * rec_field return sum( ref_class(rst_axis, dd_in=dd_in)( - sym.inverse_metric_derivative( - rst_axis, expr.op.xyz_axis, - ambient_dim=self.ambient_dim, dim=self.dim) *rec_field) + sym.inverse_metric_derivative( + rst_axis, expr.op.xyz_axis, + ambient_dim=self.ambient_dim, dim=self.dim) * rec_field) for rst_axis in range(self.dim)) - if isinstance(expr.op, op.MassOperator): return op.RefMassOperator(expr.op.dd_in, expr.op.dd_out)( jac_in * self.rec(expr.field)) @@ -530,8 +529,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): in_inv_op = op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_in) trans_mass_op = op.RefMassOperator(expr.op.dd_in, expr.op.dd_out) out_inv_op = op.RefInverseMassOperator(expr.op.dd_out, expr.op.dd_out) - return out_inv_op(trans_mass_op( in_inv_op(self.rec(expr.field))/jac_in )) - + return out_inv_op(trans_mass_op(in_inv_op(self.rec(expr.field))/jac_in)) elif isinstance(expr.op, op.FaceMassOperator): jac_in_surf = sym.area_element(self.ambient_dim, self.dim - 1, diff --git a/test/test_grudge.py b/test/test_grudge.py index 94456bb7..97bccfbd 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -332,6 +332,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type assert eoc_rec.order_estimate() > order + @pytest.mark.parametrize("order", [3, 4, 5]) def test_convergence_curvi(ctx_factory, order): """Test whether a curvilinear mesh converges""" @@ -363,12 +364,12 @@ def test_convergence_curvi(ctx_factory, order): from grudge.models.advection import WeakAdvectionOperator analytic_sol = bind(discr, u_analytic(sym.nodes(dims))) - analytic_init = bind(discr, u_analytic(sym.nodes(dims)) ) + analytic_init = bind(discr, u_analytic(sym.nodes(dims))) fields = analytic_init(queue, t=0) op = WeakAdvectionOperator(advec_v, - inflow_u=u_analytic(sym.nodes(dims, sym.BTAG_ALL) ), + inflow_u=u_analytic(sym.nodes(dims, sym.BTAG_ALL)), flux_type="central") bound_op = bind(discr, op.sym_operator()) @@ -396,13 +397,10 @@ def test_convergence_curvi(ctx_factory, order): esc = event.state_component step += 1 - sol = analytic_sol(queue, t=step * dt) total_error = norm(queue, u=(esc - sol)) / norm(queue, u=sol) eoc_rec.add_data_point(1.0/n, total_error) - - print(eoc_rec.pretty_print(abscissa_label="h", error_label="L2 Error")) -- GitLab From 0cb91b5aeb93e0c6a7a2c49b2cc529989531d65c Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 27 Feb 2018 22:45:35 -0600 Subject: [PATCH 3/6] Fixed quadrature and added a comment --- grudge/symbolic/mappers/__init__.py | 4 +++- test/test_grudge.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index e1a26f01..1a5dd32b 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -518,7 +518,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): ref_class(rst_axis, dd_in=dd_in)( sym.inverse_metric_derivative( rst_axis, expr.op.xyz_axis, - ambient_dim=self.ambient_dim, dim=self.dim) * rec_field) + ambient_dim=self.ambient_dim, dim=self.dim, dd=dd_in) * rec_field) for rst_axis in range(self.dim)) if isinstance(expr.op, op.MassOperator): @@ -526,6 +526,8 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): jac_in * self.rec(expr.field)) elif isinstance(expr.op, op.InverseMassOperator): + # Chan, Hewett, and Warburton. "Weight-Adjusted Discontinuous + # Galerkin Methods: Curvilinear Meshes." Siam Journal on Scientific Computing in_inv_op = op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_in) trans_mass_op = op.RefMassOperator(expr.op.dd_in, expr.op.dd_out) out_inv_op = op.RefInverseMassOperator(expr.op.dd_out, expr.op.dd_out) diff --git a/test/test_grudge.py b/test/test_grudge.py index 97bccfbd..3976f246 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -334,7 +334,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type @pytest.mark.parametrize("order", [3, 4, 5]) -def test_convergence_curvi(ctx_factory, order): +def test_convergence_curvi_advection(ctx_factory, order): """Test whether a curvilinear mesh converges""" cl_ctx = cl.create_some_context() -- GitLab From ffe1fa775e39359b39f906c2c85a692d859a6d66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Wed, 28 Feb 2018 13:09:13 -0500 Subject: [PATCH 4/6] Fix article reference, add DOI link --- grudge/symbolic/mappers/__init__.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 1a5dd32b..5cfe1faf 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -526,8 +526,14 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): jac_in * self.rec(expr.field)) elif isinstance(expr.op, op.InverseMassOperator): + # follows + # # Chan, Hewett, and Warburton. "Weight-Adjusted Discontinuous - # Galerkin Methods: Curvilinear Meshes." Siam Journal on Scientific Computing + # Galerkin Methods: Curvilinear Meshes." + # SIAM Journal on Scientific Computing + # + # https://doi.org/10.1137/16M1089198 + in_inv_op = op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_in) trans_mass_op = op.RefMassOperator(expr.op.dd_in, expr.op.dd_out) out_inv_op = op.RefInverseMassOperator(expr.op.dd_out, expr.op.dd_out) -- GitLab From 433241aab3ba0e321064c35cdd8ef70930892456 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Wed, 28 Feb 2018 13:11:25 -0500 Subject: [PATCH 5/6] Flake8 fixes --- grudge/symbolic/mappers/__init__.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 5cfe1faf..9fc5c709 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -518,7 +518,8 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): ref_class(rst_axis, dd_in=dd_in)( sym.inverse_metric_derivative( rst_axis, expr.op.xyz_axis, - ambient_dim=self.ambient_dim, dim=self.dim, dd=dd_in) * rec_field) + ambient_dim=self.ambient_dim, dim=self.dim, dd=dd_in) + * rec_field) for rst_axis in range(self.dim)) if isinstance(expr.op, op.MassOperator): @@ -537,7 +538,10 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): in_inv_op = op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_in) trans_mass_op = op.RefMassOperator(expr.op.dd_in, expr.op.dd_out) out_inv_op = op.RefInverseMassOperator(expr.op.dd_out, expr.op.dd_out) - return out_inv_op(trans_mass_op(in_inv_op(self.rec(expr.field))/jac_in)) + return out_inv_op( + trans_mass_op(in_inv_op(self.rec(expr.field)) + / + jac_in)) elif isinstance(expr.op, op.FaceMassOperator): jac_in_surf = sym.area_element(self.ambient_dim, self.dim - 1, -- GitLab From 647ac75635d3e1594b18020418534902346c6318 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 6 Mar 2018 13:09:00 -0600 Subject: [PATCH 6/6] Removed whitespace in empty line to placate flake --- grudge/symbolic/mappers/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 9fc5c709..e533fcd7 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -534,7 +534,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): # SIAM Journal on Scientific Computing # # https://doi.org/10.1137/16M1089198 - + in_inv_op = op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_in) trans_mass_op = op.RefMassOperator(expr.op.dd_in, expr.op.dd_out) out_inv_op = op.RefInverseMassOperator(expr.op.dd_out, expr.op.dd_out) -- GitLab