diff --git a/examples/advection/curvi.py b/examples/advection/curvi.py new file mode 100644 index 0000000000000000000000000000000000000000..b192199e1b851a2cb5eaea047634473a19ddca74 --- /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 0000000000000000000000000000000000000000..5d73495a7a84c8f7be58f09790d1a5b8e34a6447 --- /dev/null +++ b/examples/advection/test.py @@ -0,0 +1,105 @@ +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))) + + 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) / sym.sqrt(jac_bound)), + 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 + if visualize: + final_t = final_t * 10 + + 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) + esc = J_mul(queue, u=esc) + + + 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], True) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index ca11cbf0d92c96b6c9d00a6e035b7b840c2a85cd..967e4ece3d43797d06ed8bedf34def7795996d92 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -508,12 +508,29 @@ 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) + 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( + 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))) #* sym.sqrt(jac_tag_o) + return sum( sym.inverse_metric_derivative( rst_axis, expr.op.xyz_axis, @@ -526,14 +543,19 @@ 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)( - 1/jac_in * self.rec(expr.field)) + self.rec(expr.field)) # * sym.sqrt(jac_tag) 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)# * 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)) elif isinstance(expr.op, op.StiffnessOperator): return op.RefMassOperator()(jac_noquad * @@ -548,7 +570,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(