From c426bbd8698ca3feb78cc9c37d52f370d6491741 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sat, 20 Jan 2018 16:00:34 -0600 Subject: [PATCH 1/3] Got a curvi test-like thing working --- examples/advection/curvi.py | 118 ++++++++++++++++++++++++++++++++++++ examples/advection/test.py | 95 +++++++++++++++++++++++++++++ 2 files changed, 213 insertions(+) create mode 100644 examples/advection/curvi.py create mode 100644 examples/advection/test.py diff --git a/examples/advection/curvi.py b/examples/advection/curvi.py new file mode 100644 index 00000000..b192199e --- /dev/null +++ b/examples/advection/curvi.py @@ -0,0 +1,118 @@ +# Copyright (C) 2007 Andreas Kloeckner +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + + +import numpy as np +import pyopencl as cl # noqa +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest # noqa + +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + + + +def main(write_output=True, order=4): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + + from meshmode.mesh.generation import generate_warped_rect_mesh + mesh = generate_warped_rect_mesh(2, order, 20) + + + dt_factor = 4 + h = 1/20 + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + + c = np.array([0.1,0.1]) + norm_c = la.norm(c) + + + flux_type = "central" + + + def f(x): + return sym.sin(3*x) + + def u_analytic(x): + return f(-np.dot(c, x)/norm_c+sym.var("t", sym.DD_SCALAR)*norm_c) + + from grudge.models.advection import WeakAdvectionOperator + from meshmode.mesh import BTAG_ALL, BTAG_NONE + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + op = WeakAdvectionOperator(c, + inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)), + flux_type=flux_type) + + bound_op = bind(discr, op.sym_operator()) + + u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0) + + def rhs(t, u): + return bound_op(queue, t=t, u=u) + + final_time = 0.3 + + dt = dt_factor * h/order**2 + nsteps = (final_time // dt) + 1 + dt = final_time/nsteps + 1e-15 + + + from grudge.shortcuts import set_up_rk4 + dt_stepper = set_up_rk4("u", dt, u, rhs) + + last_u = None + + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + + step = 0 + + norm = bind(discr, sym.norm(2, sym.var("u"))) + + for event in dt_stepper.run(t_end=final_time): + if isinstance(event, dt_stepper.StateComputed): + + step += 1 + + #print(step, event.t, norm(queue, u=event.state_component[0])) + vis.write_vtk_file("fld-%04d.vtu" % step, + [ ("u", event.state_component) ]) + + + + + + + +if __name__ == "__main__": + main() + + diff --git a/examples/advection/test.py b/examples/advection/test.py new file mode 100644 index 00000000..8fab9eb7 --- /dev/null +++ b/examples/advection/test.py @@ -0,0 +1,95 @@ +import numpy as np +import pyopencl as cl # noqa + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +import numpy.linalg as la + + +def test_convergence_maxwell(order, ns, visualize=False): + """Test whether 3D maxwells actually converges""" + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + dims = 2 + for n in ns: + h = 1.0/n + dt_factor = 4 + + #from meshmode.mesh.generation import generate_regular_rect_mesh + #mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + #n=(n, n), order=order) + + 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))) + fields = analytic_sol(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 + + 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) + + print("dt=%g nsteps=%d" % (dt, nsteps)) + + 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 + print(step) + + 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) + + if visualize: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + vis.write_vtk_file("fld-thing.vtu", [ ("a", esc), ("s", sol), ("e", (esc-sol)) ]) + return + + + + print(eoc_rec.pretty_print(abscissa_label="h", + error_label="L2 Error")) + + #assert eoc_rec.order_estimate() > order + + +if __name__ == "__main__": + test_convergence_maxwell(3, [15,20]) -- GitLab From 4f6c4e2cd7eae467286cfab88fd3aef6568d5afe Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sun, 11 Feb 2018 19:47:21 -0600 Subject: [PATCH 2/3] WIP: Somehow, curvi does not immediately explode --- examples/advection/test.py | 4 +++- grudge/symbolic/mappers/__init__.py | 28 ++++++++++++++++++++++++---- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/examples/advection/test.py b/examples/advection/test.py index 8fab9eb7..734385f9 100644 --- a/examples/advection/test.py +++ b/examples/advection/test.py @@ -52,6 +52,8 @@ def test_convergence_maxwell(order, ns, visualize=False): return bound_op(queue, t=t, u=u) final_t = 0.1 + if visualize: + final_t = final_t * 10 dt = dt_factor * h/order**2 nsteps = (final_t // dt) + 1 @@ -92,4 +94,4 @@ def test_convergence_maxwell(order, ns, visualize=False): if __name__ == "__main__": - test_convergence_maxwell(3, [15,20]) + test_convergence_maxwell(3, [15,20], True) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index ca11cbf0..2b47e33b 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -508,12 +508,28 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): jac_noquad = sym.area_element(self.ambient_dim, dim, dd=expr.op.dd_in.with_qtag(sym.QTAG_NONE)) - def rewrite_derivative(ref_class, field, dd_in, with_jacobian=True): + def rewrite_derivative(ref_class, field, dd_in, with_jacobian=True, SUPERBAD_HACK_FLAG = False): jac_tag = sym.area_element(self.ambient_dim, self.dim, dd=dd_in) rec_field = self.rec(field) if with_jacobian: + just_rec_field = rec_field # SUPERBAD rec_field = jac_tag * rec_field + if SUPERBAD_HACK_FLAG: + 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) * + just_rec_field) + for rst_axis in range(self.dim)) - \ + sum(sym.RefMassOperator(expr.op.dd_in, expr.op.dd_out)( + sym.inverse_metric_derivative( + rst_axis, expr.op.xyz_axis, + ambient_dim=self.ambient_dim, dim=self.dim) * + just_rec_field * sym.RefDiffOperator(rst_axis)(jac_tag) / (2 * jac_tag)) + for rst_axis in range(self.dim)) + return sum( sym.inverse_metric_derivative( rst_axis, expr.op.xyz_axis, @@ -527,13 +543,17 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): elif isinstance(expr.op, op.InverseMassOperator): return op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_out)( - 1/jac_in * self.rec(expr.field)) + self.rec(expr.field)) elif isinstance(expr.op, op.FaceMassOperator): + from grudge.symbolic.primitives import DOFDesc + jac_vol = sym.interp(DOFDesc("vol"), expr.op.dd_in)(sym.area_element(self.ambient_dim, self.dim, dd=DOFDesc("vol"))) jac_in_surf = sym.area_element(self.ambient_dim, self.dim - 1, dd=expr.op.dd_in) return op.RefFaceMassOperator(expr.op.dd_in, expr.op.dd_out)( - jac_in_surf * self.rec(expr.field)) + jac_in_surf * self.rec(expr.field)/ jac_vol) + #return op.RefFaceMassOperator(expr.op.dd_in, expr.op.dd_out)( + #jac_in_surf * self.rec(expr.field)) elif isinstance(expr.op, op.StiffnessOperator): return op.RefMassOperator()(jac_noquad * @@ -548,7 +568,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): elif isinstance(expr.op, op.StiffnessTOperator): return rewrite_derivative( op.RefStiffnessTOperator, - expr.field, dd_in=expr.op.dd_in) + expr.field, dd_in=expr.op.dd_in, with_jacobian=True, SUPERBAD_HACK_FLAG=True) elif isinstance(expr.op, op.MInvSTOperator): return self.rec( -- GitLab From c8aa09278d5f461e91cdae53c78eb63308d955ff Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 12 Feb 2018 17:56:45 -0600 Subject: [PATCH 3/3] WIP: Attempt to compensate with Jacobians --- examples/advection/test.py | 12 ++++++++++-- grudge/symbolic/mappers/__init__.py | 10 ++++++---- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/examples/advection/test.py b/examples/advection/test.py index 734385f9..5d73495a 100644 --- a/examples/advection/test.py +++ b/examples/advection/test.py @@ -40,10 +40,16 @@ def test_convergence_maxwell(order, ns, visualize=False): from grudge.models.advection import WeakAdvectionOperator analytic_sol = bind(discr, u_analytic(sym.nodes(dims))) - fields = analytic_sol(queue, t=0) + jac_tag = sym.area_element(2, 2, dd=sym.DOFDesc("vol")) + analytic_init = bind(discr, u_analytic(sym.nodes(dims)) / sym.sqrt(jac_tag)) + J_mul = bind(discr, sym.var("u", sym.DOFDesc("vol")) * sym.sqrt(jac_tag)) + + fields = analytic_init(queue, t=0) + + jac_bound = sym.interp(sym.DOFDesc("vol"), sym.BTAG_ALL)(jac_tag) op = WeakAdvectionOperator(advec_v, - inflow_u=u_analytic(sym.nodes(dims, sym.BTAG_ALL)), + inflow_u=u_analytic(sym.nodes(dims, sym.BTAG_ALL) / sym.sqrt(jac_bound)), flux_type="central") bound_op = bind(discr, op.sym_operator()) @@ -74,6 +80,8 @@ def test_convergence_maxwell(order, ns, visualize=False): step += 1 print(step) + esc = J_mul(queue, u=esc) + sol = analytic_sol(queue, t=step * dt) total_error = norm(queue, u=(esc - sol)) / norm(queue, u=sol) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 2b47e33b..967e4ece 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -510,13 +510,14 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): def rewrite_derivative(ref_class, field, dd_in, with_jacobian=True, SUPERBAD_HACK_FLAG = False): jac_tag = sym.area_element(self.ambient_dim, self.dim, dd=dd_in) + jac_tag_o = sym.area_element(self.ambient_dim, self.dim, dd=expr.op.dd_out) rec_field = self.rec(field) if with_jacobian: just_rec_field = rec_field # SUPERBAD rec_field = jac_tag * rec_field if SUPERBAD_HACK_FLAG: - return sum( + return (sum( ref_class(rst_axis, dd_in=dd_in)( sym.inverse_metric_derivative( rst_axis, expr.op.xyz_axis, @@ -528,7 +529,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): rst_axis, expr.op.xyz_axis, ambient_dim=self.ambient_dim, dim=self.dim) * just_rec_field * sym.RefDiffOperator(rst_axis)(jac_tag) / (2 * jac_tag)) - for rst_axis in range(self.dim)) + for rst_axis in range(self.dim))) #* sym.sqrt(jac_tag_o) return sum( sym.inverse_metric_derivative( @@ -542,8 +543,9 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): jac_in * self.rec(expr.field)) elif isinstance(expr.op, op.InverseMassOperator): + jac_tag = sym.area_element(self.ambient_dim, self.dim, dd=expr.op.dd_out) return op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_out)( - self.rec(expr.field)) + self.rec(expr.field)) # * sym.sqrt(jac_tag) elif isinstance(expr.op, op.FaceMassOperator): from grudge.symbolic.primitives import DOFDesc @@ -551,7 +553,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): jac_in_surf = sym.area_element(self.ambient_dim, self.dim - 1, dd=expr.op.dd_in) return op.RefFaceMassOperator(expr.op.dd_in, expr.op.dd_out)( - jac_in_surf * self.rec(expr.field)/ jac_vol) + jac_in_surf * self.rec(expr.field)/ jac_vol)# * sym.sqrt(sym.area_element(self.ambient_dim, self.dim, dd=DOFDesc("vol"))) #return op.RefFaceMassOperator(expr.op.dd_in, expr.op.dd_out)( #jac_in_surf * self.rec(expr.field)) -- GitLab