From c14939c3e9d22b8e31044f649ff8d848d0359e7f Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sat, 24 Feb 2018 02:22:00 -0600 Subject: [PATCH 01/15] 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 02/15] 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 03/15] 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 fcf4a86c7b40f5ea635bcdbe9953562cd06ed66e Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 21 Mar 2018 15:17:17 -0500 Subject: [PATCH 04/15] Everything is broken --- examples/gas_dynamics/euler/sine-wave.py | 230 +++++++------------ examples/gas_dynamics/euler/test.py | 106 +++++++++ grudge/models/gas_dynamics/__init__.py | 272 +++++++++++++---------- grudge/symbolic/mappers/__init__.py | 63 +++++- grudge/symbolic/operators.py | 4 + grudge/symbolic/primitives.py | 1 + 6 files changed, 412 insertions(+), 264 deletions(-) create mode 100644 examples/gas_dynamics/euler/test.py diff --git a/examples/gas_dynamics/euler/sine-wave.py b/examples/gas_dynamics/euler/sine-wave.py index 451befdb..576722b1 100644 --- a/examples/gas_dynamics/euler/sine-wave.py +++ b/examples/gas_dynamics/euler/sine-wave.py @@ -13,17 +13,10 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +import pyopencl as cl # noqa +import numpy as np - - -from __future__ import division -from __future__ import absolute_import -from __future__ import print_function -import numpy -import numpy.linalg as la - - - +from grudge import sym, bind, DGDiscretizationWithBoundaries class SineWave: def __init__(self): @@ -32,155 +25,102 @@ class SineWave: self.prandtl = 0.72 self.spec_gas_const = 287.1 - def __call__(self, t, x_vec): - rho = 2 + numpy.sin(x_vec[0] + x_vec[1] + x_vec[2] - 2 * t) - velocity = numpy.array([1, 1, 0]) + def sym_operator(self): + t = sym.var("t", sym.DD_SCALAR) + sym_nds = sym.nodes(3) + rho = 2 + sym.sin(sym_nds[0] + sym_nds[1] + sym_nds[2] - 2 * t) + velocity = np.array([1, 1, 0]) p = 1 - e = p/(self.gamma-1) + rho/2 * numpy.dot(velocity, velocity) + e = p/(self.gamma-1) + rho/2 * np.dot(velocity, velocity) rho_u = rho * velocity[0] rho_v = rho * velocity[1] rho_w = rho * velocity[2] - from grudge.tools import join_fields + from pytools.obj_array import join_fields return join_fields(rho, e, rho_u, rho_v, rho_w) def properties(self): return(self.gamma, self.mu, self.prandtl, self.spec_gas_const) - def volume_interpolant(self, t, discr): - return discr.convert_volume( - self(t, discr.nodes.T), - kind=discr.compute_kind) - - def boundary_interpolant(self, t, discr, tag): - return discr.convert_boundary( - self(t, discr.get_boundary(tag).nodes.T), - tag=tag, kind=discr.compute_kind) - - - +def main(order = 3, visualize=False): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) -def main(final_time=1, write_output=False): - from grudge.backends import guess_run_context - rcon = guess_run_context() - from grudge.tools import EOCRecorder, to_obj_array + from pytools.convergence import EOCRecorder eoc_rec = EOCRecorder() - if rcon.is_head_rank: - from grudge.mesh import make_box_mesh - mesh = make_box_mesh((0,0,0), (10,10,10), max_volume=0.5) - mesh_data = rcon.distribute_mesh(mesh) - else: - mesh_data = rcon.receive_mesh() - - for order in [3, 4, 5]: - discr = rcon.make_discretization(mesh_data, order=order, - default_scalar_type=numpy.float64) - - from grudge.visualization import SiloVisualizer, VtkVisualizer - vis = VtkVisualizer(discr, rcon, "sinewave-%d" % order) - #vis = SiloVisualizer(discr, rcon) - - sinewave = SineWave() - fields = sinewave.volume_interpolant(0, discr) - gamma, mu, prandtl, spec_gas_const = sinewave.properties() - - from grudge.mesh import BTAG_ALL - from grudge.models.gas_dynamics import GasDynamicsOperator - op = GasDynamicsOperator(dimensions=mesh.dimensions, gamma=gamma, mu=mu, - prandtl=prandtl, spec_gas_const=spec_gas_const, - bc_inflow=sinewave, bc_outflow=sinewave, bc_noslip=sinewave, - inflow_tag=BTAG_ALL, source=None) - - euler_ex = op.bind(discr) - - max_eigval = [0] - def rhs(t, q): - ode_rhs, speed = euler_ex(t, q) - max_eigval[0] = speed - return ode_rhs - rhs(0, fields) - - if rcon.is_head_rank: - print("---------------------------------------------") - print("order %d" % order) - print("---------------------------------------------") - print("#elements=", len(mesh.elements)) - - from grudge.timestep import RK4TimeStepper - stepper = RK4TimeStepper() - - # diagnostics setup --------------------------------------------------- - from pytools.log import LogManager, add_general_quantities, \ - add_simulation_quantities, add_run_info - - if write_output: - log_name = ("euler-sinewave-%(order)d-%(els)d.dat" - % {"order":order, "els":len(mesh.elements)}) - else: - log_name = False - logmgr = LogManager(log_name, "w", rcon.communicator) - add_run_info(logmgr) - add_general_quantities(logmgr) - add_simulation_quantities(logmgr) - discr.add_instrumentation(logmgr) - stepper.add_instrumentation(logmgr) - - logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"]) - - # timestep loop ------------------------------------------------------- - try: - from grudge.timestep import times_and_steps - step_it = times_and_steps( - final_time=final_time, logmgr=logmgr, - max_dt_getter=lambda t: op.estimate_timestep(discr, - stepper=stepper, t=t, max_eigenvalue=max_eigval[0])) - - for step, t, dt in step_it: - #if step % 10 == 0: - if write_output: - visf = vis.make_file("sinewave-%d-%04d" % (order, step)) - - #from pyvisfile.silo import DB_VARTYPE_VECTOR - vis.add_data(visf, - [ - ("rho", discr.convert_volume(op.rho(fields), kind="numpy")), - ("e", discr.convert_volume(op.e(fields), kind="numpy")), - ("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")), - ("u", discr.convert_volume(op.u(fields), kind="numpy")), - - #("true_rho", op.rho(true_fields)), - #("true_e", op.e(true_fields)), - #("true_rho_u", op.rho_u(true_fields)), - #("true_u", op.u(true_fields)), - - #("rhs_rho", op.rho(rhs_fields)), - #("rhs_e", op.e(rhs_fields)), - #("rhs_rho_u", op.rho_u(rhs_fields)), - ], - #expressions=[ - #("diff_rho", "rho-true_rho"), - #("diff_e", "e-true_e"), - #("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR), - - #("p", "0.4*(e- 0.5*(rho_u*u))"), - #], - time=t, step=step - ) - visf.close() - - fields = stepper(fields, t, dt, rhs) - - finally: - vis.close() - logmgr.close() - discr.close() - - true_fields = sinewave.volume_interpolant(t, discr) - eoc_rec.add_data_point(order, discr.norm(fields-true_fields)) - print() - print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) + n = 6 + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5, -0.5), b=(0.5, 0.5, 0.5), + n=(n, n, n), order=order) + + from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory # noqa + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, + quad_tag_to_group_factory={ + #"product": None, + "product": QuadratureSimplexGroupFactory(order=3*order) + }) + + + sinewave = SineWave() + bound_sine = bind(discr, sinewave.sym_operator()) + fields = bound_sine(queue, t=0) + gamma, mu, prandtl, spec_gas_const = sinewave.properties() + + from meshmode.mesh import BTAG_ALL + from grudge.models.gas_dynamics import GasDynamicsOperator + op = GasDynamicsOperator(dimensions=3, gamma=gamma, mu=mu, + prandtl=prandtl, spec_gas_const=spec_gas_const, + bc_inflow=sinewave, bc_outflow=sinewave, bc_noslip=sinewave, + inflow_tag=BTAG_ALL, source=None) + + euler_ex = bind(discr, op.sym_operator()) + + max_eigval = [0] + def rhs(t, q): + ode_rhs, speed = euler_ex(t, q) + max_eigval[0] = speed + return ode_rhs + rhs(0, fields) + + from grudge.shortcuts import set_up_rk4 + stepper = set_up_rk4("q", dt, fields, rhs) + + if visualize: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + + step = 0 + for event in dt_stepper.run(t_end=final_t): + if isinstance(event, dt_stepper.StateComputed): + assert event.component_id == "u" + fields = event.state_component + + vis.write_vtk_file("sinewave-%d-%04d" % (order, step), + [ + ("rho", op.rho(fields)), + ("e", op.e(fields)), + ("rho_u", op.rho_u(fields)), + ("u", op.u(fields)), + + #("true_rho", op.rho(true_fields)), + #("true_e", op.e(true_fields)), + #("true_rho_u", op.rho_u(true_fields)), + #("true_u", op.u(true_fields)), + + #("rhs_rho", op.rho(rhs_fields)), + #("rhs_e", op.e(rhs_fields)), + #("rhs_rho_u", op.rho_u(rhs_fields)), + ]) + + step += 1 + print(step) + + true_fields = bound_sine(queue, t=t) + eoc_rec.add_data_point(order, discr.norm(fields-true_fields)) + print() + print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) diff --git a/examples/gas_dynamics/euler/test.py b/examples/gas_dynamics/euler/test.py new file mode 100644 index 00000000..44d478f1 --- /dev/null +++ b/examples/gas_dynamics/euler/test.py @@ -0,0 +1,106 @@ +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/models/gas_dynamics/__init__.py b/grudge/models/gas_dynamics/__init__.py index e5a8ddc9..d76a8fd2 100644 --- a/grudge/models/gas_dynamics/__init__.py +++ b/grudge/models/gas_dynamics/__init__.py @@ -32,36 +32,21 @@ THE SOFTWARE. import numpy import grudge.tools -import grudge.mesh -import grudge.data +import meshmode.mesh from grudge.models import TimeDependentOperator from pytools import Record from grudge.tools import is_zero -from grudge.second_order import ( - StabilizedCentralSecondDerivative, - CentralSecondDerivative, - IPDGSecondDerivative) -from grudge.symbolic.primitives import make_common_subexpression as cse +from pymbolic.primitives import make_common_subexpression as cse from pytools import memoize_method -from grudge.symbolic.tools import make_sym_vector from pytools.obj_array import make_obj_array, join_fields +from grudge import sym AXES = ["x", "y", "z", "w"] -from grudge.symbolic.operators import ( - QuadratureGridUpsampler, - QuadratureInteriorFacesGridUpsampler) -to_vol_quad = QuadratureGridUpsampler("gasdyn_vol") - -# It is recommended (though not required) that these two -# remain the same so that they can be computed together -# by the CUDA backend - -to_int_face_quad = QuadratureInteriorFacesGridUpsampler("gasdyn_face") -to_bdry_quad = QuadratureGridUpsampler("gasdyn_face") +# TODO noslip? # {{{ equations of state @@ -141,8 +126,9 @@ class GasDynamicsOperator(TimeDependentOperator): supersonic_inflow_tag="supersonic_inflow", supersonic_outflow_tag="supersonic_outflow", source=None, - second_order_scheme=CentralSecondDerivative(), + second_order_scheme=None, artificial_viscosity_mode=None, + quad_tag = "product" ): """ :param source: should implement @@ -151,9 +137,6 @@ class GasDynamicsOperator(TimeDependentOperator): :param artificial_viscosity_mode: """ - from grudge.data import ( - TimeConstantGivenFunction, - ConstantGivenFunction) if gamma is not None: if equation_of_state is not None: @@ -165,9 +148,7 @@ class GasDynamicsOperator(TimeDependentOperator): equation_of_state = GammaLawEOS(gamma) - dull_bc = TimeConstantGivenFunction( - ConstantGivenFunction(make_obj_array( - [1, 1] + [0]*dimensions))) + dull_bc = make_obj_array( [1, 1] + [0]*dimensions) if bc_inflow is None: bc_inflow = dull_bc if bc_outflow is None: @@ -205,6 +186,8 @@ class GasDynamicsOperator(TimeDependentOperator): raise ValueError("artificial_viscosity_mode has an invalid value") self.artificial_viscosity_mode = artificial_viscosity_mode + self.quad_tag = quad_tag + @@ -212,14 +195,19 @@ class GasDynamicsOperator(TimeDependentOperator): # {{{ conversions --------------------------------------------------------- def state(self): - return make_sym_vector("q", self.dimensions+2) + return sym.make_sym_array("q", self.dimensions+2) @memoize_method def volq_state(self): + quad_dd = sym.DOFDesc("vol", self.quad_tag) + to_vol_quad = sym.interp("vol", quad_dd) return cse(to_vol_quad(self.state()), "vol_quad_state") @memoize_method def faceq_state(self): + int_faces_dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR, self.quad_tag) + to_int_face_quad = sym.interp("vol", int_faces_dd) + #return cse(sym.int_tpair(self.state(), self.quad_tag), "face_quad_state") return cse(to_int_face_quad(self.state()), "face_quad_state") @memoize_method @@ -298,7 +286,7 @@ class GasDynamicsOperator(TimeDependentOperator): def primitive_to_conservative(self, prims, use_cses=True): if not use_cses: - from grudge.symbolic.primitives import make_common_subexpression as cse + def cse(x, name): return sym.cse(x,name) else: def cse(x, name): return x @@ -314,7 +302,7 @@ class GasDynamicsOperator(TimeDependentOperator): def conservative_to_primitive(self, q, use_cses=True): if use_cses: - from grudge.symbolic.primitives import make_common_subexpression as cse + def cse(x, name): return sym.cse(x,name) else: def cse(x, name): return x @@ -334,10 +322,10 @@ class GasDynamicsOperator(TimeDependentOperator): "sound_speed") u = self.cse_u(state) speed = cse(sqrt(numpy.dot(u, u)), "norm_u") + sound_speed - return ElementwiseMaxOperator()(speed) + return ElementwiseMaxOperator("vol","vol")(speed) def bind_characteristic_velocity(self, discr): - state = make_sym_vector("q", self.dimensions+2) + state = sym.make_sym_array("q", self.dimensions+2) compiled = discr.compile( self.characteristic_velocity_optemplate(state)) @@ -360,10 +348,10 @@ class GasDynamicsOperator(TimeDependentOperator): int_flux_operand=faceq_var, bdry_flux_int_operand=faceq_var) - self.second_order_scheme.grad(grad_tgt, - bc_getter=self.get_boundary_condition_for, - dirichlet_tags=self.get_boundary_tags(), - neumann_tags=[]) + #self.second_order_scheme.grad(grad_tgt, + #bc_getter=self.get_boundary_condition_for, + #dirichlet_tags=self.get_boundary_tags(), + #neumann_tags=[]) return grad_tgt.minv_all @@ -391,7 +379,7 @@ class GasDynamicsOperator(TimeDependentOperator): # {{{ viscous stress tensor def tau(self, to_quad_op, state, mu=None): - faceq_state = self.faceq_state() + #faceq_state = self.faceq_state() dimensions = self.dimensions @@ -488,8 +476,10 @@ class GasDynamicsOperator(TimeDependentOperator): BC is linearized. """ if state0 is None: - state0 = make_sym_vector(bc_name, self.dimensions+2) + state0 = sym.make_sym_array(bc_name, self.dimensions+2) + boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) + to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd) state0 = cse(to_bdry_quad(state0)) rho0 = self.rho(state0) @@ -498,8 +488,7 @@ class GasDynamicsOperator(TimeDependentOperator): c0 = (self.equation_of_state.gamma * p0 / rho0)**0.5 - from grudge.symbolic import RestrictToBoundary - bdrize_op = RestrictToBoundary(tag) + bdrize_op = sym.interp("vol",sym.BTAG_ALL) class SingleBCInfo(Record): pass @@ -516,8 +505,8 @@ class GasDynamicsOperator(TimeDependentOperator): - p0, "dpm")) def outflow_state(self, state): - from grudge.symbolic import make_normal - normal = make_normal(self.outflow_tag, self.dimensions) + + normal = sym.normal(sym.DTAG_BOUNDARY(self.outflow_tag), self.dimensions) bc = self.make_bc_info("bc_q_out", self.outflow_tag, state) # see grudge/doc/maxima/euler.mac @@ -552,29 +541,25 @@ class GasDynamicsOperator(TimeDependentOperator): + normal*numpy.dot(normal, bc.dumvec)/2 + bc.dpm*normal/(2*bc.c0*bc.rho0), "bc_u_"+name)) def inflow_state(self, state): - from grudge.symbolic import make_normal - normal = make_normal(self.inflow_tag, self.dimensions) + normal = sym.normal(sym.DTAG_BOUNDARY(self.inflow_tag), self.dimensions) bc = self.make_bc_info("bc_q_in", self.inflow_tag, state) return self.inflow_state_inner(normal, bc, "inflow") def noslip_state(self, state): - from grudge.symbolic import make_normal state0 = join_fields( - make_sym_vector("bc_q_noslip", 2), + sym.make_sym_array("bc_q_noslip", 2), [0]*self.dimensions) - normal = make_normal(self.noslip_tag, self.dimensions) + normal = sym.normal(sym.DTAG_BOUNDARY(self.noslip_tag), self.dimensions) bc = self.make_bc_info("bc_q_noslip", self.noslip_tag, state, state0) return self.inflow_state_inner(normal, bc, "noslip") def wall_state(self, state): - from grudge.symbolic import RestrictToBoundary - bc = RestrictToBoundary(self.wall_tag)(state) + bc = sym.interp("vol",sym.BTAG_ALL)(state) wall_rho = self.rho(bc) wall_e = self.e(bc) # <3 eve wall_rho_u = self.rho_u(bc) - from grudge.symbolic import make_normal - normal = make_normal(self.wall_tag, self.dimensions) + normal = sym.normal(sym.DTAG_BOUNDARY(self.wall_tag), self.dimensions) return join_fields( wall_rho, @@ -596,13 +581,12 @@ class GasDynamicsOperator(TimeDependentOperator): def get_conservative_boundary_conditions(self): state = self.state() - from grudge.symbolic import RestrictToBoundary return { self.supersonic_inflow_tag: - make_sym_vector("bc_q_supersonic_in", self.dimensions+2), + sym.make_sym_array("bc_q_supersonic_in", self.dimensions+2), self.supersonic_outflow_tag: - RestrictToBoundary(self.supersonic_outflow_tag)( - (state)), + #sym.interp("vol",self.supersonic_outflow_tag)(state), + sym.interp("vol",sym.BTAG_ALL)(state), self.wall_tag: self.wall_state(state), } @@ -615,11 +599,10 @@ class GasDynamicsOperator(TimeDependentOperator): def _normalize_expr(self, expr): """Normalize expressions for use as hash keys.""" from grudge.symbolic.mappers import ( - QuadratureUpsamplerRemover, + QuadratureRemover, CSERemover) - return CSERemover()( - QuadratureUpsamplerRemover({}, do_warn=False)(expr)) + return CSERemover()(QuadratureRemover()(expr)) @memoize_method def _get_norm_primitive_exprs(self): @@ -632,6 +615,8 @@ class GasDynamicsOperator(TimeDependentOperator): def get_boundary_condition_for(self, tag, expr): prim_bcs = self.get_primitive_boundary_conditions() cons_bcs = self.get_conservative_boundary_conditions() + boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) + to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd) if tag in prim_bcs: # BC is given in primitive variables, avoid converting @@ -652,23 +637,21 @@ class GasDynamicsOperator(TimeDependentOperator): # 'cbstate' is the boundary state in conservative variables. - from grudge.symbolic.mappers import QuadratureUpsamplerRemover - expr = QuadratureUpsamplerRemover({}, do_warn=False)(expr) + from grudge.symbolic.mappers import QuadratureRemover + expr = QuadratureRemover()(expr) def subst_func(expr): from pymbolic.primitives import Subscript, Variable + boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) + to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd) if isinstance(expr, Subscript): assert (isinstance(expr.aggregate, Variable) and expr.aggregate.name == "q") return cbstate[expr.index] - elif isinstance(expr, Variable) and expr.name =="sensor": - from grudge.symbolic import RestrictToBoundary - result = RestrictToBoundary(tag)(self.sensor()) - return cse(to_bdry_quad(result), "bdry_sensor") - from grudge.symbolic import SubstitutionMapper + from grudge.symbolic.mappers import SubstitutionMapper return SubstitutionMapper(subst_func)(expr) # }}} @@ -681,35 +664,37 @@ class GasDynamicsOperator(TimeDependentOperator): operand=vol_operand, int_flux_operand=int_face_operand) - self.second_order_scheme.div(div_tgt, - bc_getter=self.get_boundary_condition_for, - dirichlet_tags=list(self.get_boundary_tags()), - neumann_tags=[]) + #self.second_order_scheme.div(div_tgt, + #bc_getter=self.get_boundary_condition_for, + #dirichlet_tags=list(self.get_boundary_tags()), + #neumann_tags=[]) return div_tgt.minv_all - def make_second_order_part(self): - state = self.state() - faceq_state = self.faceq_state() - volq_state = self.volq_state() - - volq_tau_mat = self.tau(to_vol_quad, state) - faceq_tau_mat = self.tau(to_int_face_quad, state) - - return join_fields( - 0, - self.div( - numpy.sum(volq_tau_mat*self.cse_u(volq_state), axis=1) - + self.heat_conduction_grad(to_vol_quad) - , - numpy.sum(faceq_tau_mat*self.cse_u(faceq_state), axis=1) - + self.heat_conduction_grad(to_int_face_quad) - , - ), - [ - self.div(volq_tau_mat[i], faceq_tau_mat[i]) - for i in range(self.dimensions)] - ) + #def make_second_order_part(self): + #state = self.state() + #faceq_state = self.faceq_state() + #volq_state = self.volq_state() + #int_faces_dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR, self.quad_tag) + #to_int_face_quad = sym.interp("vol", int_faces_dd) +# + #volq_tau_mat = self.tau(to_vol_quad, state) + #faceq_tau_mat = self.tau(to_int_face_quad, state) +# + #return join_fields( + #0, + #self.div( + #numpy.sum(volq_tau_mat*self.cse_u(volq_state), axis=1) + #+ self.heat_conduction_grad(to_vol_quad) + #, + #numpy.sum(faceq_tau_mat*self.cse_u(faceq_state), axis=1) + #+ self.heat_conduction_grad(to_int_face_quad) + #, + #), + #[ + #self.div(volq_tau_mat[i], faceq_tau_mat[i]) + #for i in range(self.dimensions)] + #) # }}} @@ -717,25 +702,77 @@ class GasDynamicsOperator(TimeDependentOperator): def make_extra_terms(self): return 0 + def make_lax_friedrichs_flux(self, wave_speed, state, fluxes, bdry_tags_states_and_fluxes,strong): + #from hedge.flux import FluxVectorPlaceholder, flux_max + + def lff(fvph): + n = len(state) + d = len(fluxes) + normal = sym.normal(fvph.dd, d) + #fvph = FluxVectorPlaceholder(len(state)*(1+d)+1) + + wave_speed_ph = fvph[0] + state_ph = fvph[1:1+n] + fluxes_ph = [fvph[1+i*n:1+(i+1)*n] for i in range(1, d+1)] + + from pymbolic.primitives import Max + penalty = Max((wave_speed_ph.int,wave_speed_ph.ext))*(state_ph.ext-state_ph.int) + + if not strong: + num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) + - penalty) + else: + num_flux = 0.5*(sum(n_i*(f_i.int-f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) + + penalty) + return num_flux + + all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) + def flux(pair): + return sym.interp(pair.dd, all_faces_dd)(lff(pair)) + + #def no_int_int_pair(expression): + #from grudge.symbolic.primitives import TracePair + #i = expression + #e = cse(sym.OppositeInteriorFaceSwap()(i)) + #return TracePair(sym.DOFDesc("int_faces", self.quad_tag), i, e) + + + + #from hedge.optemplate import get_flux_operator + #flux_op = get_flux_operator(num_flux) + int_operand = join_fields(wave_speed, state, *fluxes) + + #return (flux(int_operand) + sum( + #flux(sym.bv_tpair(sym.DTAG_BOUNDARY(tag),int_operand,join_fields(0, bdry_state, *bdry_fluxes) )) + #for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) + + return (flux(sym.int_tpair(int_operand, self.quad_tag)) + sum( + flux(sym.bv_tpair(sym.DTAG_BOUNDARY(tag),int_operand,join_fields(0, bdry_state, *bdry_fluxes) )) + for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) + + + def sym_operator(self, sensor_scaling=None, viscosity_only=False): u = self.cse_u rho = self.cse_rho rho_u = self.rho_u p = self.p e = self.e + int_faces_dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR, self.quad_tag) + to_int_face_quad = sym.interp("vol", int_faces_dd) # {{{ artificial diffusion - def make_artificial_diffusion(): - if self.artificial_viscosity_mode not in ["diffusion"]: - return 0 + #def make_artificial_diffusion(): + #if self.artificial_viscosity_mode not in ["diffusion"]: + #return 0 - dq = self.grad_of_state() + #dq = self.grad_of_state() - return make_obj_array([ - self.div( - to_vol_quad(self.sensor())*to_vol_quad(dq[i]), - to_int_face_quad(self.sensor())*to_int_face_quad(dq[i])) - for i in range(dq.shape[0])]) + #return make_obj_array([ + #self.div( + #to_vol_quad(self.sensor())*to_vol_quad(dq[i]), + #to_int_face_quad(self.sensor())*to_int_face_quad(dq[i])) + #for i in range(dq.shape[0])]) # }}} # {{{ state setup @@ -751,13 +788,10 @@ class GasDynamicsOperator(TimeDependentOperator): has_viscosity = not is_zero(self.get_mu(self.state(), to_quad_op=None)) # }}} + boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) + to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd) # {{{ operator assembly ----------------------------------------------- - from grudge.flux.tools import make_lax_friedrichs_flux - from grudge.symbolic.operators import InverseMassOperator - - from grudge.symbolic.tools import make_stiffness_t - primitive_bcs_as_quad_conservative = dict( (tag, self.primitive_to_conservative(to_bdry_quad(bc))) for tag, bc in @@ -769,29 +803,31 @@ class GasDynamicsOperator(TimeDependentOperator): self.get_boundary_condition_for(tag, s_i) for s_i in state]) return tag, bc, self.flux(bc) - first_order_part = InverseMassOperator()( - numpy.dot(make_stiffness_t(self.dimensions), volq_flux) - - make_lax_friedrichs_flux( - wave_speed=cse(to_int_face_quad(speed), "emax_c"), - state=self.faceq_state(), fluxes=faceq_flux, - bdry_tags_states_and_fluxes=[ - get_bc_tuple(tag) for tag in self.get_boundary_tags()], - strong=False)) + + first_order_part = sym.InverseMassOperator()( + numpy.dot(sym.stiffness_t(self.dimensions), volq_flux) + - self.make_lax_friedrichs_flux( + cse(to_int_face_quad(speed), "emax_c"), + self.faceq_state(), + faceq_flux, + [get_bc_tuple(tag) for tag in self.get_boundary_tags()], + False) + ) if viscosity_only: first_order_part = 0*first_order_part result = join_fields( - first_order_part - + self.make_second_order_part() - + make_artificial_diffusion() - + self.make_extra_terms(), + first_order_part , + #+ self.make_second_order_part() + #+ make_artificial_diffusion() + #+ self.make_extra_terms(), speed) if self.source is not None: result = result + join_fields( - make_sym_vector("source_vect", len(self.state())), + sym.make_sym_array("source_vect", len(self.state())), # extra field for speed 0) @@ -812,7 +848,7 @@ class GasDynamicsOperator(TimeDependentOperator): sensor_scaling=sensor_scaling, viscosity_only=False)) - from grudge.mesh import check_bc_coverage + from meshmode.mesh import check_bc_coverage check_bc_coverage(discr.mesh, [ self.inflow_tag, self.outflow_tag, diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 1a5dd32b..1a254b4d 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -183,6 +183,7 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin): map_nodal_sum = map_elementwise_linear map_nodal_min = map_elementwise_linear map_nodal_max = map_elementwise_linear + map_elementwise_max = map_elementwise_linear map_stiffness = map_elementwise_linear map_diff = map_elementwise_linear @@ -722,7 +723,13 @@ class NoCSEStringifyMapper(StringifyMapper): # {{{ quadrature support class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper): - """Checks whether all quadratu + """ Checks that all quadrature tags exist in qttg (quad_tag_to_group_factory), + and then removes all quadrature changes that are None in qttg. + + For example: if you have an expression that interpolates to a QTAG called + "product", and the discretization is made with qttg = {"product": None}, + then all operators will be changed to ignore "product" QTAG. + """ def __init__(self, quad_tag_to_group_factory): @@ -733,6 +740,7 @@ class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper): map_common_subexpression_uncached = \ IdentityMapper.map_common_subexpression + # def _process_dd(self, dd, location_descr): from grudge.symbolic.primitives import DOFDesc, QTAG_NONE if dd.quadrature_tag is not QTAG_NONE: @@ -783,6 +791,58 @@ class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper): dd = self._process_dd(expr.dd, location_descr=type(expr).__name__) return type(expr)(expr.axis, dd) +class QuadratureRemover(CSECachingMapperMixin, IdentityMapper): + """ Removes all quadrature use + """ + + def __init__(self): + IdentityMapper.__init__(self) + CSECachingMapperMixin.__init__(self) + + map_common_subexpression_uncached = \ + IdentityMapper.map_common_subexpression + + def _process_dd(self, dd): + from grudge.symbolic.primitives import DOFDesc, QTAG_NONE + return DOFDesc(dd.domain_tag, QTAG_NONE) + + def map_operator_binding(self, expr): + dd_in_orig = dd_in = expr.op.dd_in + dd_out_orig = dd_out = expr.op.dd_out + + dd_in = self._process_dd(dd_in) + dd_out = self._process_dd(dd_out) + + if dd_in_orig == dd_in and dd_out_orig == dd_out: + # unchanged + return IdentityMapper.map_operator_binding(self, expr) + + import grudge.symbolic.operators as op + # changed + + if dd_in == dd_out and isinstance(expr.op, op.InterpolationOperator): + # This used to be to-quad interpolation and has become a no-op. + # Remove it. + return self.rec(expr.field) + + if isinstance(expr.op, op.StiffnessTOperator): + new_op = type(expr.op)(expr.op.xyz_axis, dd_in, dd_out) + elif isinstance(expr.op, (op.FaceMassOperator, op.InterpolationOperator)): + new_op = type(expr.op)(dd_in, dd_out) + else: + raise NotImplementedError("do not know how to modify dd_in and dd_out " + "in %s" % type(expr.op).__name__) + + return new_op(self.rec(expr.field)) + + #def map_ones(self, expr): + #dd = self._process_dd(expr.dd, location_descr=type(expr).__name__) + #return type(expr)(dd) + + #def map_node_coordinate_component(self, expr): + #dd = self._process_dd(expr.dd, location_descr=type(expr).__name__) + #return type(expr)(expr.axis, dd) + # }}} @@ -1083,6 +1143,7 @@ class CollectorMixin(OperatorReducerMixin, LocalOpReducerMixin, FluxOpReducerMix map_node_coordinate_component = map_grudge_variable map_operator = map_grudge_variable + map_elementwise_max = map_grudge_variable #TODO: lolwhat why? # I'm not sure this works still. diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 194b0716..e5404bb6 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -57,6 +57,10 @@ class Operator(pymbolic.primitives.Expression): def __init__(self, dd_in, dd_out): self.dd_in = _sym().as_dofdesc(dd_in) self.dd_out = _sym().as_dofdesc(dd_out) + #TODO: remove this. According to Python docs this creates weird circular references + # And a bad time for everyone + import inspect + self.stack = inspect.stack() def stringifier(self): from grudge.symbolic.mappers import StringifyMapper diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 00f7b8d3..f61931ec 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -311,6 +311,7 @@ class Variable(HasDOFDesc, ExpressionBase, pymbolic.primitives.Variable): HasDOFDesc.__init__(self, dd) pymbolic.primitives.Variable.__init__(self, name) + def __getinitargs__(self): return (self.name, self.dd,) -- GitLab From 38c06cb895ded43294e54ac170fc8825002aa9a1 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 24 Apr 2018 18:32:20 -0500 Subject: [PATCH 05/15] euler attempt --- examples/gas_dynamics/euler/sine-wave.py | 2 +- examples/wave/wave-min.py | 8 ++-- grudge/models/gas_dynamics/__init__.py | 59 +++++++++++++++--------- grudge/symbolic/dofdesc_inference.py | 4 ++ grudge/symbolic/mappers/__init__.py | 5 +- grudge/symbolic/operators.py | 2 +- grudge/symbolic/primitives.py | 7 +++ 7 files changed, 59 insertions(+), 28 deletions(-) diff --git a/examples/gas_dynamics/euler/sine-wave.py b/examples/gas_dynamics/euler/sine-wave.py index 576722b1..d6e49445 100644 --- a/examples/gas_dynamics/euler/sine-wave.py +++ b/examples/gas_dynamics/euler/sine-wave.py @@ -79,7 +79,7 @@ def main(order = 3, visualize=False): max_eigval = [0] def rhs(t, q): - ode_rhs, speed = euler_ex(t, q) + ode_rhs, speed = euler_ex(queue, t=t, q=q) max_eigval[0] = speed return ode_rhs rhs(0, fields) diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py index 89f4bfa8..e366a992 100644 --- a/examples/wave/wave-min.py +++ b/examples/wave/wave-min.py @@ -35,12 +35,12 @@ def main(write_output=True, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - dims = 3 + dims = 2 from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh( a=(-0.5,)*dims, b=(0.5,)*dims, - n=(16,)*dims) + n=(9,)*dims) if mesh.dim == 2: dt = 0.04 @@ -114,7 +114,7 @@ def main(write_output=True, order=4): print(step, event.t, norm(queue, u=event.state_component[0]), time()-t_last_step) - if step % 10 == 0: + if step % 3 == 0: vis.write_vtk_file("fld-%04d.vtu" % step, [ ("u", event.state_component[0]), @@ -124,4 +124,4 @@ def main(write_output=True, order=4): if __name__ == "__main__": - main() + main(order = 4) diff --git a/grudge/models/gas_dynamics/__init__.py b/grudge/models/gas_dynamics/__init__.py index d76a8fd2..e1bfc191 100644 --- a/grudge/models/gas_dynamics/__init__.py +++ b/grudge/models/gas_dynamics/__init__.py @@ -46,7 +46,6 @@ AXES = ["x", "y", "z", "w"] -# TODO noslip? # {{{ equations of state @@ -164,6 +163,8 @@ class GasDynamicsOperator(TimeDependentOperator): self.spec_gas_const = spec_gas_const self.mu = mu + self.BAD_DEBUG = None + self.bc_inflow = bc_inflow self.bc_outflow = bc_outflow self.bc_noslip = bc_noslip @@ -479,8 +480,11 @@ class GasDynamicsOperator(TimeDependentOperator): state0 = sym.make_sym_array(bc_name, self.dimensions+2) boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) - to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd) - state0 = cse(to_bdry_quad(state0)) + boundary_dd = sym.DOFDesc(sym.DTAG_BOUNDARY(tag), self.quad_tag) + to_bdry_quad = sym.interp(sym.BTAG_ALL, boundary_dd) + #state0 = cse(to_bdry_quad(state0)) + to_bdry_quad_f_a = sym.interp("vol", boundary_dd) + state0 = cse(to_bdry_quad_f_a(state0)) rho0 = self.rho(state0) p0 = self.cse_p(state0) @@ -488,7 +492,8 @@ class GasDynamicsOperator(TimeDependentOperator): c0 = (self.equation_of_state.gamma * p0 / rho0)**0.5 - bdrize_op = sym.interp("vol",sym.BTAG_ALL) + #bdrize_op = sym.interp("vol",sym.BTAG_ALL) + bdrize_op = sym.interp("vol",boundary_dd) class SingleBCInfo(Record): pass @@ -585,8 +590,9 @@ class GasDynamicsOperator(TimeDependentOperator): self.supersonic_inflow_tag: sym.make_sym_array("bc_q_supersonic_in", self.dimensions+2), self.supersonic_outflow_tag: - #sym.interp("vol",self.supersonic_outflow_tag)(state), - sym.interp("vol",sym.BTAG_ALL)(state), + #sym.interp("vol",self.supersonic_outflow_tag)(state), #TODO: It should be this according to my em translation <- but (it spits out an error) + sym.interp("vol",sym.DTAG_BOUNDARY(self.supersonic_outflow_tag))(state), #TODO: Attempt 2 + #sym.interp("vol",sym.BTAG_ALL)(state), self.wall_tag: self.wall_state(state), } @@ -633,7 +639,8 @@ class GasDynamicsOperator(TimeDependentOperator): # BC is given in conservative variables, no potential # for optimization. - cbstate = to_bdry_quad(cons_bcs[tag]) + #cbstate = to_bdry_quad(cons_bcs[tag]) + cbstate = cons_bcs[tag] # 'cbstate' is the boundary state in conservative variables. @@ -642,8 +649,8 @@ class GasDynamicsOperator(TimeDependentOperator): def subst_func(expr): from pymbolic.primitives import Subscript, Variable - boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) - to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd) + #boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) + #to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd) if isinstance(expr, Subscript): assert (isinstance(expr.aggregate, Variable) @@ -708,26 +715,27 @@ class GasDynamicsOperator(TimeDependentOperator): def lff(fvph): n = len(state) d = len(fluxes) - normal = sym.normal(fvph.dd, d) + normal = sym.normal(fvph.dd, d) #TODO: Also contradiction between this and some dd errors at the next todo #fvph = FluxVectorPlaceholder(len(state)*(1+d)+1) wave_speed_ph = fvph[0] state_ph = fvph[1:1+n] fluxes_ph = [fvph[1+i*n:1+(i+1)*n] for i in range(1, d+1)] + print(fvph.dd) from pymbolic.primitives import Max penalty = Max((wave_speed_ph.int,wave_speed_ph.ext))*(state_ph.ext-state_ph.int) if not strong: - num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) + num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) # TODO dd error between flux and norm here - penalty) else: - num_flux = 0.5*(sum(n_i*(f_i.int-f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) + num_flux = 0.5*(sum(n_i*(f_i.int-f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) # TODO (or here) + penalty) return num_flux all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) - def flux(pair): + def fluxy(pair): return sym.interp(pair.dd, all_faces_dd)(lff(pair)) #def no_int_int_pair(expression): @@ -746,8 +754,14 @@ class GasDynamicsOperator(TimeDependentOperator): #flux(sym.bv_tpair(sym.DTAG_BOUNDARY(tag),int_operand,join_fields(0, bdry_state, *bdry_fluxes) )) #for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) - return (flux(sym.int_tpair(int_operand, self.quad_tag)) + sum( - flux(sym.bv_tpair(sym.DTAG_BOUNDARY(tag),int_operand,join_fields(0, bdry_state, *bdry_fluxes) )) + + + #return (flux(sym.int_tpair(int_operand, self.quad_tag)) + sum( + #flux(sym.bv_tpair(sym.DTAG_BOUNDARY(tag),int_operand,join_fields(0, bdry_state, *bdry_fluxes) )) + #for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) + + return (fluxy(sym.int_tpair(int_operand, self.quad_tag)) + sum( + fluxy(sym.bv_tpair(sym.DOFDesc(sym.DTAG_BOUNDARY(tag), self.quad_tag),int_operand,join_fields(0, bdry_state, *bdry_fluxes) )) for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) @@ -804,15 +818,18 @@ class GasDynamicsOperator(TimeDependentOperator): return tag, bc, self.flux(bc) + quad_dd = sym.DOFDesc("vol", self.quad_tag) + all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) + stiff_t = sym.stiffness_t(self.dimensions, quad_dd, "vol") first_order_part = sym.InverseMassOperator()( - numpy.dot(sym.stiffness_t(self.dimensions), volq_flux) - - self.make_lax_friedrichs_flux( - cse(to_int_face_quad(speed), "emax_c"), - self.faceq_state(), - faceq_flux, + numpy.dot(stiff_t, volq_flux) + - sym.FaceMassOperator(all_faces_dd, "vol")(self.make_lax_friedrichs_flux( + speed, + self.state(), + self.flux(self.state()), #faceq_flux, [get_bc_tuple(tag) for tag in self.get_boundary_tags()], - False) + False)) ) if viscosity_only: diff --git a/grudge/symbolic/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py index 832c6a03..379969bc 100644 --- a/grudge/symbolic/dofdesc_inference.py +++ b/grudge/symbolic/dofdesc_inference.py @@ -157,6 +157,10 @@ class DOFDescInferenceMapper(RecursiveMapper, CSECachingMapperMixin): def map_nodal_sum(self, expr, enclosing_prec): return DOFDesc(DTAG_SCALAR) + def map_max(self, expr): + return self.map_multi_child(expr, expr.children) + + map_nodal_max = map_nodal_sum map_nodal_min = map_nodal_sum diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 1a254b4d..59e4b084 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -135,6 +135,7 @@ class OperatorReducerMixin(LocalOpReducerMixin, FluxOpReducerMixin): map_nodal_sum = _map_op_base map_nodal_min = _map_op_base map_nodal_max = _map_op_base + map_elementwise_max = _map_op_base map_stiffness = _map_op_base map_diff = _map_op_base @@ -615,6 +616,9 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper): def map_nodal_min(self, expr, enclosing_prec): return "NodalMin" + self._format_op_dd(expr) + def map_elementwise_max(self, expr, enclosing_prec): + return "ElementwiseMax" + self._format_op_dd(expr) + # }}} # {{{ global differentiation @@ -1143,7 +1147,6 @@ class CollectorMixin(OperatorReducerMixin, LocalOpReducerMixin, FluxOpReducerMix map_node_coordinate_component = map_grudge_variable map_operator = map_grudge_variable - map_elementwise_max = map_grudge_variable #TODO: lolwhat why? # I'm not sure this works still. diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index e5404bb6..fcf87d23 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -117,7 +117,7 @@ interp = InterpolationOperator class ElementwiseMaxOperator(Operator): mapper_method = intern("map_elementwise_max") - + # {{{ nodal reduction: sum, integral, max diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index f61931ec..9fd7b790 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -198,6 +198,13 @@ class DOFDesc(object): self.domain_tag = domain_tag self.quadrature_tag = quadrature_tag + + #TODO: remove this. According to Python docs this creates weird circular references + # And a bad time for everyone + import inspect + self.stack = inspect.stack() + + def is_scalar(self): return self.domain_tag is DTAG_SCALAR -- GitLab From 62e728629a761fc101e59425b6576e363d74ff35 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 16 May 2018 00:07:14 -0500 Subject: [PATCH 06/15] euler schenanigans --- grudge/execution.py | 31 ++++++++++++++++++++++++++----- grudge/symbolic/compiler.py | 5 +++-- 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/grudge/execution.py b/grudge/execution.py index 2da2e508..f8b69459 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -233,15 +233,36 @@ class ExecutionMapper(mappers.Evaluator, return result def map_elementwise_max(self, op, field_expr): - from grudge._internal import perform_elwise_max - field = self.rec(field_expr) + pu.db + inp = self.rec(field_expr) + import pymbolic.primitives as p + var = p.Variable + + i = var("i") + @memoize_in(self.bound_op, "face_elwisemax_knl") + def knl(): + knl = lp.make_kernel( + "{[i,j]: 0<=i Date: Wed, 16 May 2018 23:42:56 -0500 Subject: [PATCH 07/15] Removed stack saves --- grudge/symbolic/compiler.py | 4 ++++ grudge/symbolic/operators.py | 4 ++-- grudge/symbolic/primitives.py | 4 ++-- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index fb6ba40c..163d2d2b 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -494,6 +494,10 @@ class Code(object): if pre_assign_check is not None: pre_assign_check(target, value) + #if isinstance(value, Pyopenclarray): + #v2 = value.get(queue) + #if np.isnan(v2).any(): + #kaboom context[target] = value futures.extend(new_futures) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index fcf87d23..31dc3166 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -59,8 +59,8 @@ class Operator(pymbolic.primitives.Expression): self.dd_out = _sym().as_dofdesc(dd_out) #TODO: remove this. According to Python docs this creates weird circular references # And a bad time for everyone - import inspect - self.stack = inspect.stack() + #import inspect + #self.stack = inspect.stack() def stringifier(self): from grudge.symbolic.mappers import StringifyMapper diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 9fd7b790..a29e481a 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -201,8 +201,8 @@ class DOFDesc(object): #TODO: remove this. According to Python docs this creates weird circular references # And a bad time for everyone - import inspect - self.stack = inspect.stack() + #import inspect + #self.stack = inspect.stack() def is_scalar(self): -- GitLab From e23a8dcb65dc87965e228f3aebd88e1f30a325fe Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Sun, 20 May 2018 17:32:46 -0500 Subject: [PATCH 08/15] more debugging --- examples/gas_dynamics/euler/sine-wave.py | 47 +++++++++++--------- grudge/execution.py | 55 ++++++++++++++++++++---- grudge/models/gas_dynamics/__init__.py | 8 ++-- grudge/symbolic/compiler.py | 10 +++-- 4 files changed, 85 insertions(+), 35 deletions(-) diff --git a/examples/gas_dynamics/euler/sine-wave.py b/examples/gas_dynamics/euler/sine-wave.py index d6e49445..d596d81f 100644 --- a/examples/gas_dynamics/euler/sine-wave.py +++ b/examples/gas_dynamics/euler/sine-wave.py @@ -74,18 +74,23 @@ def main(order = 3, visualize=False): prandtl=prandtl, spec_gas_const=spec_gas_const, bc_inflow=sinewave, bc_outflow=sinewave, bc_noslip=sinewave, inflow_tag=BTAG_ALL, source=None) + print(sym.pretty(op.sym_operator())) euler_ex = bind(discr, op.sym_operator()) + print(euler_ex) max_eigval = [0] def rhs(t, q): - ode_rhs, speed = euler_ex(queue, t=t, q=q) - max_eigval[0] = speed - return ode_rhs + #ode_rhs, speed = euler_ex(queue, t=t, q=q) + things = euler_ex(queue, t=t, q=q) + max_eigval[0] = things[-1] + return things[0:-1] rhs(0, fields) from grudge.shortcuts import set_up_rk4 - stepper = set_up_rk4("q", dt, fields, rhs) + dt = 0.000001 + final_t = 0.0001 + dt_stepper = set_up_rk4("q", dt, fields, rhs) if visualize: from grudge.shortcuts import make_visualizer @@ -94,25 +99,27 @@ def main(order = 3, visualize=False): step = 0 for event in dt_stepper.run(t_end=final_t): if isinstance(event, dt_stepper.StateComputed): - assert event.component_id == "u" + #assert event.component_id == "u" + print(event.component_id) fields = event.state_component - vis.write_vtk_file("sinewave-%d-%04d" % (order, step), - [ - ("rho", op.rho(fields)), - ("e", op.e(fields)), - ("rho_u", op.rho_u(fields)), - ("u", op.u(fields)), + if visualize: + vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, step), + [ + ("rho", op.rho(fields)), + ("e", op.e(fields)), + ("rho_u", op.rho_u(fields)), + ("u", op.u(fields)), - #("true_rho", op.rho(true_fields)), - #("true_e", op.e(true_fields)), - #("true_rho_u", op.rho_u(true_fields)), - #("true_u", op.u(true_fields)), + #("true_rho", op.rho(true_fields)), + #("true_e", op.e(true_fields)), + #("true_rho_u", op.rho_u(true_fields)), + #("true_u", op.u(true_fields)), - #("rhs_rho", op.rho(rhs_fields)), - #("rhs_e", op.e(rhs_fields)), - #("rhs_rho_u", op.rho_u(rhs_fields)), - ]) + #("rhs_rho", op.rho(rhs_fields)), + #("rhs_e", op.e(rhs_fields)), + #("rhs_rho_u", op.rho_u(rhs_fields)), + ]) step += 1 print(step) @@ -126,7 +133,7 @@ def main(order = 3, visualize=False): if __name__ == "__main__": - main() + main(visualize=True) diff --git a/grudge/execution.py b/grudge/execution.py index f8b69459..451f483b 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -114,6 +114,7 @@ class ExecutionMapper(mappers.Evaluator, cached_name += "abs" else: cached_name += expr.function.name + print(cached_name + ":") @memoize_in(self.bound_op, cached_name) def knl(): @@ -128,6 +129,8 @@ class ExecutionMapper(mappers.Evaluator, assert len(args) == 1 evt, (out,) = knl()(self.queue, a=args[0]) + print(out) + return out def map_nodal_sum(self, op, field_expr): @@ -233,8 +236,41 @@ class ExecutionMapper(mappers.Evaluator, return result def map_elementwise_max(self, op, field_expr): - pu.db + field = self.rec(field_expr) + + from grudge.tools import is_zero + if is_zero(field): + return 0 + + @memoize_in(self.bound_op, "elwise_max_knl") + def knl(): + knl = lp.make_kernel( + """{[k,i,j]: + 0<=k Date: Sun, 20 May 2018 23:00:31 -0500 Subject: [PATCH 09/15] added prints for kernel calls --- grudge/execution.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/grudge/execution.py b/grudge/execution.py index 451f483b..c5af8fa1 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -391,6 +391,10 @@ class ExecutionMapper(mappers.Evaluator, kwargs["grdg_n"] = discr.nnodes evt, result_dict = kdescr.loopy_kernel(self.queue, **kwargs) + if "expr_217" in result_dict.keys(): + knl = lp.set_options(kdescr.loopy_kernel, "write_cl") + knl(self.queue, **kwargs) + pu.db return list(result_dict.items()), [] def map_insn_assign(self, insn): -- GitLab From 3bf487fab6fb651703f08a52482fa00df2ddb2ff Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 28 May 2018 01:40:41 -0500 Subject: [PATCH 10/15] added vortex.py support --- examples/gas_dynamics/euler/sine-wave.py | 4 +- examples/gas_dynamics/euler/vortex.py | 191 ++++++------------ .../gas_dynamics/gas_dynamics_initials.py | 28 ++- grudge/execution.py | 3 +- grudge/symbolic/operators.py | 4 +- grudge/symbolic/primitives.py | 4 +- 6 files changed, 87 insertions(+), 147 deletions(-) diff --git a/examples/gas_dynamics/euler/sine-wave.py b/examples/gas_dynamics/euler/sine-wave.py index d596d81f..13aeb587 100644 --- a/examples/gas_dynamics/euler/sine-wave.py +++ b/examples/gas_dynamics/euler/sine-wave.py @@ -50,7 +50,7 @@ def main(order = 3, visualize=False): from pytools.convergence import EOCRecorder eoc_rec = EOCRecorder() - n = 6 + n = 3 from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh(a=(-0.5, -0.5, -0.5), b=(0.5, 0.5, 0.5), n=(n, n, n), order=order) @@ -133,7 +133,7 @@ def main(order = 3, visualize=False): if __name__ == "__main__": - main(visualize=True) + main(order=2,visualize=True) diff --git a/examples/gas_dynamics/euler/vortex.py b/examples/gas_dynamics/euler/vortex.py index 9f9703f4..2bf9f8a6 100644 --- a/examples/gas_dynamics/euler/vortex.py +++ b/examples/gas_dynamics/euler/vortex.py @@ -21,62 +21,55 @@ from __future__ import absolute_import from __future__ import print_function import numpy +import pyopencl as cl # noqa +from grudge import sym, bind, DGDiscretizationWithBoundaries - -def main(write_output=True): +def main(write_output=True, visualize=False): from pytools import add_python_path_relative_to_script add_python_path_relative_to_script("..") + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) - from grudge.backends import guess_run_context - rcon = guess_run_context() - from grudge.tools import EOCRecorder + from pytools.convergence import EOCRecorder eoc_rec = EOCRecorder() - if rcon.is_head_rank: - from grudge.mesh.generator import \ - make_rect_mesh, \ - make_centered_regular_rect_mesh + from meshmode.mesh.generation import generate_regular_rect_mesh - refine = 4 - mesh = make_centered_regular_rect_mesh((0,-5), (10,5), n=(9,9), - post_refine_factor=refine) - mesh_data = rcon.distribute_mesh(mesh) - else: - mesh_data = rcon.receive_mesh() for order in [3, 4, 5]: + refine = 4 + mesh = generate_regular_rect_mesh(a=(0,-5), b=(10,5), n=(4,4), order=order) from gas_dynamics_initials import Vortex flow = Vortex() from grudge.models.gas_dynamics import ( GasDynamicsOperator, PolytropeEOS, GammaLawEOS) - from grudge.mesh import BTAG_ALL + from meshmode.mesh import BTAG_ALL + from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory # works equally well for GammaLawEOS op = GasDynamicsOperator(dimensions=2, mu=flow.mu, prandtl=flow.prandtl, spec_gas_const=flow.spec_gas_const, equation_of_state=PolytropeEOS(flow.gamma), - bc_inflow=flow, bc_outflow=flow, bc_noslip=flow, + bc_inflow=flow.volume_interpolant(), bc_outflow=flow.volume_interpolant(), bc_noslip=flow.volume_interpolant(), inflow_tag=BTAG_ALL, source=None) - discr = rcon.make_discretization(mesh_data, order=order, - default_scalar_type=numpy.float64, - quad_min_degrees={ - "gasdyn_vol": 3*order, - "gasdyn_face": 3*order, - }, - tune_for=op.sym_operator(), - debug=["cuda_no_plan"]) + discr = DGDiscretizationWithBoundaries(cl_ctx,mesh, order=order, + quad_tag_to_group_factory={ + "product": QuadratureSimplexGroupFactory(order=3*order) + }) - from grudge.visualization import SiloVisualizer, VtkVisualizer - vis = VtkVisualizer(discr, rcon, "vortex-%d" % order) - #vis = SiloVisualizer(discr, rcon) + if visualize: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) - fields = flow.volume_interpolant(0, discr) - euler_ex = op.bind(discr) + bound_flow_vol = bind(discr, flow.volume_interpolant()) + fields = bound_flow_vol(queue, t=0) + + euler_ex = bind(discr, op.sym_operator()) max_eigval = [0] def rhs(t, q): @@ -85,92 +78,53 @@ def main(write_output=True): return ode_rhs rhs(0, fields) - if rcon.is_head_rank: - print("---------------------------------------------") - print("order %d" % order) - print("---------------------------------------------") - print("#elements=", len(mesh.elements)) - # limiter ------------------------------------------------------------ - from grudge.models.gas_dynamics import SlopeLimiter1NEuler - limiter = SlopeLimiter1NEuler(discr, flow.gamma, 2, op) + #from grudge.models.gas_dynamics import SlopeLimiter1NEuler + #limiter = SlopeLimiter1NEuler(discr, flow.gamma, 2, op) - from grudge.timestep.runge_kutta import SSP3TimeStepper + #from grudge.timestep.runge_kutta import SSP3TimeStepper #stepper = SSP3TimeStepper(limiter=limiter) - stepper = SSP3TimeStepper( - vector_primitive_factory=discr.get_vector_primitive_factory()) - - #from grudge.timestep import RK4TimeStepper - #stepper = RK4TimeStepper() - - # diagnostics setup --------------------------------------------------- - from pytools.log import LogManager, add_general_quantities, \ - add_simulation_quantities, add_run_info - - if write_output: - log_file_name = "euler-%d.dat" % order - else: - log_file_name = None - - logmgr = LogManager(log_file_name, "w", rcon.communicator) - add_run_info(logmgr) - add_general_quantities(logmgr) - add_simulation_quantities(logmgr) - discr.add_instrumentation(logmgr) - stepper.add_instrumentation(logmgr) - - logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"]) - - # timestep loop ------------------------------------------------------- - try: - final_time = flow.final_time - from grudge.timestep import times_and_steps - step_it = times_and_steps( - final_time=final_time, logmgr=logmgr, - max_dt_getter=lambda t: op.estimate_timestep(discr, - stepper=stepper, t=t, max_eigenvalue=max_eigval[0])) - - print("run until t=%g" % final_time) - for step, t, dt in step_it: - if step % 10 == 0 and write_output: - #if False: - visf = vis.make_file("vortex-%d-%04d" % (order, step)) - - #true_fields = vortex.volume_interpolant(t, discr) - - from pyvisfile.silo import DB_VARTYPE_VECTOR - vis.add_data(visf, - [ - ("rho", discr.convert_volume(op.rho(fields), kind="numpy")), - ("e", discr.convert_volume(op.e(fields), kind="numpy")), - ("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")), - ("u", discr.convert_volume(op.u(fields), kind="numpy")), - - #("true_rho", discr.convert_volume(op.rho(true_fields), kind="numpy")), - #("true_e", discr.convert_volume(op.e(true_fields), kind="numpy")), - #("true_rho_u", discr.convert_volume(op.rho_u(true_fields), kind="numpy")), - #("true_u", discr.convert_volume(op.u(true_fields), kind="numpy")), - - #("rhs_rho", discr.convert_volume(op.rho(rhs_fields), kind="numpy")), - #("rhs_e", discr.convert_volume(op.e(rhs_fields), kind="numpy")), - #("rhs_rho_u", discr.convert_volume(op.rho_u(rhs_fields), kind="numpy")), - ], - #expressions=[ - #("diff_rho", "rho-true_rho"), - #("diff_e", "e-true_e"), - #("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR), - - #("p", "0.4*(e- 0.5*(rho_u*u))"), - #], - time=t, step=step - ) - visf.close() - - fields = stepper(fields, t, dt, rhs) - #fields = limiter(fields) - - assert not numpy.isnan(numpy.sum(fields[0])) + #stepper = SSP3TimeStepper( + #vector_primitive_factory=discr.get_vector_primitive_factory()) + from grudge.shortcuts import set_up_rk4 + dt = 0.000001 + final_t = 0.0001 + dt_stepper = set_up_rk4("q", dt, fields, rhs) + + if visualize: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + + step = 0 + for event in dt_stepper.run(t_end=final_t): + if isinstance(event, dt_stepper.StateComputed): + #assert event.component_id == "u" + print(event.component_id) + fields = event.state_component + + if visualize: + vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, step), + [ + ("rho", op.rho(fields)), + ("e", op.e(fields)), + ("rho_u", op.rho_u(fields)), + ("u", op.u(fields)), + + #("true_rho", op.rho(true_fields)), + #("true_e", op.e(true_fields)), + #("true_rho_u", op.rho_u(true_fields)), + #("true_u", op.u(true_fields)), + + #("rhs_rho", op.rho(rhs_fields)), + #("rhs_e", op.e(rhs_fields)), + #("rhs_rho_u", op.rho_u(rhs_fields)), + ]) + step += 1 + print(step) + + + pu.db true_fields = flow.volume_interpolant(final_time, discr) l2_error = discr.norm(fields-true_fields) @@ -180,22 +134,9 @@ def main(write_output=True): l2_error_u = discr.norm(op.u(fields)-op.u(true_fields)) eoc_rec.add_data_point(order, l2_error) - print() print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) - logmgr.set_constant("l2_error", l2_error) - logmgr.set_constant("l2_error_rho", l2_error_rho) - logmgr.set_constant("l2_error_e", l2_error_e) - logmgr.set_constant("l2_error_rhou", l2_error_rhou) - logmgr.set_constant("l2_error_u", l2_error_u) - logmgr.set_constant("refinement", refine) - - finally: - if write_output: - vis.close() - logmgr.close() - discr.close() # after order loop assert eoc_rec.estimate_order_of_convergence()[0,1] > 6 diff --git a/examples/gas_dynamics/gas_dynamics_initials.py b/examples/gas_dynamics/gas_dynamics_initials.py index 05ab91a3..dd26417a 100644 --- a/examples/gas_dynamics/gas_dynamics_initials.py +++ b/examples/gas_dynamics/gas_dynamics_initials.py @@ -21,6 +21,7 @@ from __future__ import absolute_import import numpy import numpy.linalg as la from six.moves import range +from grudge import sym @@ -131,7 +132,10 @@ class Vortex: self.prandtl = 0.72 self.spec_gas_const = 287.1 - def __call__(self, t, x_vec): + def __call__(self): + x_vec = sym.nodes(2) + t = sym.var("t", sym.DD_SCALAR) + vortex_loc = self.center + t*self.velocity # coordinates relative to vortex center @@ -154,11 +158,6 @@ class Vortex: from grudge.tools import join_fields return join_fields(rho, e, rho*u, rho*v) - def volume_interpolant(self, t, discr): - return discr.convert_volume( - self(t, discr.nodes.T - .astype(discr.default_scalar_type)), - kind=discr.compute_kind) def boundary_interpolant(self, t, discr, tag): return discr.convert_boundary( @@ -184,7 +183,9 @@ class Vortex: self.prandtl = 0.72 self.spec_gas_const = 287.1 - def __call__(self, t, x_vec): + def __call__(self): + x_vec = sym.nodes(2) + t = sym.var("t", sym.DD_SCALAR) vortex_loc = self.center + t*self.velocity # coordinates relative to vortex center @@ -195,8 +196,8 @@ class Vortex: # also JSH/TW Nodal DG Methods, p. 209 from math import pi - r = numpy.sqrt(x_rel**2+y_rel**2) - expterm = self.beta*numpy.exp(1-r**2) + r = sym.sqrt(x_rel**2+y_rel**2) + expterm = self.beta*sym.exp(1-r**2) u = self.velocity[0] - expterm*y_rel/(2*pi) v = self.velocity[1] + expterm*x_rel/(2*pi) rho = (1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1)) @@ -204,14 +205,11 @@ class Vortex: e = p/(self.gamma-1) + rho/2*(u**2+v**2) - from grudge.tools import join_fields + from pytools.obj_array import join_fields return join_fields(rho, e, rho*u, rho*v) - def volume_interpolant(self, t, discr): - return discr.convert_volume( - self(t, discr.nodes.T - .astype(discr.default_scalar_type)), - kind=discr.compute_kind) + def volume_interpolant(self): + return self() def boundary_interpolant(self, t, discr, tag): return discr.convert_boundary( diff --git a/grudge/execution.py b/grudge/execution.py index c5af8fa1..81603b5c 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -390,6 +390,7 @@ class ExecutionMapper(mappers.Evaluator, % (name, type(v))) kwargs["grdg_n"] = discr.nnodes + kdescr.loopy_kernel = lp.set_options(kdescr.loopy_kernel, trace_assignment_values=True) evt, result_dict = kdescr.loopy_kernel(self.queue, **kwargs) if "expr_217" in result_dict.keys(): knl = lp.set_options(kdescr.loopy_kernel, "write_cl") @@ -552,7 +553,7 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, sym_operator = mappers.EmptyFluxKiller(discrwb.mesh)(sym_operator) dumper("before-cfold", sym_operator) - sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator) + #sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator) dumper("before-qcheck", sym_operator) sym_operator = mappers.QuadratureCheckerAndRemover( diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 31dc3166..fcf87d23 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -59,8 +59,8 @@ class Operator(pymbolic.primitives.Expression): self.dd_out = _sym().as_dofdesc(dd_out) #TODO: remove this. According to Python docs this creates weird circular references # And a bad time for everyone - #import inspect - #self.stack = inspect.stack() + import inspect + self.stack = inspect.stack() def stringifier(self): from grudge.symbolic.mappers import StringifyMapper diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index a29e481a..9fd7b790 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -201,8 +201,8 @@ class DOFDesc(object): #TODO: remove this. According to Python docs this creates weird circular references # And a bad time for everyone - #import inspect - #self.stack = inspect.stack() + import inspect + self.stack = inspect.stack() def is_scalar(self): -- GitLab From 08d87600714edd3004dedaa66be62f13d3e3c51f Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 30 May 2018 00:16:12 -0500 Subject: [PATCH 11/15] Mistakes were made and this contains some heat equation attempts --- examples/gas_dynamics/euler/vortex.py | 95 +++++++------- .../gas_dynamics/gas_dynamics_initials.py | 2 +- examples/heat/heat2.py | 124 ++++++++++++++++++ grudge/execution.py | 2 +- grudge/models/diffusion.py | 17 ++- grudge/models/gas_dynamics/__init__.py | 2 +- grudge/models/poisson.py | 29 ++-- grudge/models/second_order.py | 66 +++++----- 8 files changed, 233 insertions(+), 104 deletions(-) create mode 100644 examples/heat/heat2.py diff --git a/examples/gas_dynamics/euler/vortex.py b/examples/gas_dynamics/euler/vortex.py index 2bf9f8a6..7131f1b5 100644 --- a/examples/gas_dynamics/euler/vortex.py +++ b/examples/gas_dynamics/euler/vortex.py @@ -68,6 +68,7 @@ def main(write_output=True, visualize=False): bound_flow_vol = bind(discr, flow.volume_interpolant()) fields = bound_flow_vol(queue, t=0) + print(fields) euler_ex = bind(discr, op.sym_operator()) @@ -88,53 +89,53 @@ def main(write_output=True, visualize=False): #stepper = SSP3TimeStepper( #vector_primitive_factory=discr.get_vector_primitive_factory()) from grudge.shortcuts import set_up_rk4 - dt = 0.000001 - final_t = 0.0001 - dt_stepper = set_up_rk4("q", dt, fields, rhs) - - if visualize: - from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr, vis_order=order) - - step = 0 - for event in dt_stepper.run(t_end=final_t): - if isinstance(event, dt_stepper.StateComputed): - #assert event.component_id == "u" - print(event.component_id) - fields = event.state_component - - if visualize: - vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, step), - [ - ("rho", op.rho(fields)), - ("e", op.e(fields)), - ("rho_u", op.rho_u(fields)), - ("u", op.u(fields)), - - #("true_rho", op.rho(true_fields)), - #("true_e", op.e(true_fields)), - #("true_rho_u", op.rho_u(true_fields)), - #("true_u", op.u(true_fields)), - - #("rhs_rho", op.rho(rhs_fields)), - #("rhs_e", op.e(rhs_fields)), - #("rhs_rho_u", op.rho_u(rhs_fields)), - ]) - step += 1 - print(step) - - - pu.db - - true_fields = flow.volume_interpolant(final_time, discr) - l2_error = discr.norm(fields-true_fields) - l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields)) - l2_error_e = discr.norm(op.e(fields)-op.e(true_fields)) - l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields)) - l2_error_u = discr.norm(op.u(fields)-op.u(true_fields)) - - eoc_rec.add_data_point(order, l2_error) - print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) + dt = 0.000001 + final_t = 0.0001 + dt_stepper = set_up_rk4("q", dt, fields, rhs) + + if visualize: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + + step = 0 + for event in dt_stepper.run(t_end=final_t): + if isinstance(event, dt_stepper.StateComputed): + #assert event.component_id == "u" + print(event.component_id) + fields = event.state_component + + if visualize: + vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, step), + [ + ("rho", op.rho(fields)), + ("e", op.e(fields)), + ("rho_u", op.rho_u(fields)), + ("u", op.u(fields)), + + #("true_rho", op.rho(true_fields)), + #("true_e", op.e(true_fields)), + #("true_rho_u", op.rho_u(true_fields)), + #("true_u", op.u(true_fields)), + + #("rhs_rho", op.rho(rhs_fields)), + #("rhs_e", op.e(rhs_fields)), + #("rhs_rho_u", op.rho_u(rhs_fields)), + ]) + step += 1 + print(step) + + + pu.db + + true_fields = flow.volume_interpolant(final_time, discr) + l2_error = discr.norm(fields-true_fields) + l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields)) + l2_error_e = discr.norm(op.e(fields)-op.e(true_fields)) + l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields)) + l2_error_u = discr.norm(op.u(fields)-op.u(true_fields)) + + eoc_rec.add_data_point(order, l2_error) + print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) diff --git a/examples/gas_dynamics/gas_dynamics_initials.py b/examples/gas_dynamics/gas_dynamics_initials.py index dd26417a..de5db3d3 100644 --- a/examples/gas_dynamics/gas_dynamics_initials.py +++ b/examples/gas_dynamics/gas_dynamics_initials.py @@ -206,7 +206,7 @@ class Vortex: e = p/(self.gamma-1) + rho/2*(u**2+v**2) from pytools.obj_array import join_fields - return join_fields(rho, e, rho*u, rho*v) + return join_fields(rho+100, e+100, rho*u+100, rho*v+100) def volume_interpolant(self): return self() diff --git a/examples/heat/heat2.py b/examples/heat/heat2.py new file mode 100644 index 00000000..54306eed --- /dev/null +++ b/examples/heat/heat2.py @@ -0,0 +1,124 @@ +# 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_regular_rect_mesh + mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5), + n=(20, 20), order=order) + + 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.If( + sym.Comparison(sym.norm(x), ">", 0.1), + 1,2) + + def u_analytic(x): + return f(x) + + from grudge.models.diffusion import DiffusionOperator + from meshmode.mesh import BTAG_ALL, BTAG_NONE + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + op = DiffusionOperator(2, + dirichlet_tag="dirichlet", + dirichlet_bc=0,#bad + neumann_tag="neumann", + neumann_bc=1,#bad + ) + + + bound_op = bind(discr, op.m_sym_operator()) + + u = bind(discr, u_analytic(sym.nodes(dim)))(queue) + + 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/grudge/execution.py b/grudge/execution.py index 81603b5c..6e942f43 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -553,7 +553,7 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None, sym_operator = mappers.EmptyFluxKiller(discrwb.mesh)(sym_operator) dumper("before-cfold", sym_operator) - #sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator) + sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator) dumper("before-qcheck", sym_operator) sym_operator = mappers.QuadratureCheckerAndRemover( diff --git a/grudge/models/diffusion.py b/grudge/models/diffusion.py index 04c5cf57..81319a92 100644 --- a/grudge/models/diffusion.py +++ b/grudge/models/diffusion.py @@ -30,18 +30,17 @@ THE SOFTWARE. import numpy -import grudge.data from grudge.models import TimeDependentOperator from grudge.models.poisson import LaplacianOperatorBase -from grudge.second_order import CentralSecondDerivative +from grudge.models.second_order import CentralSecondDerivative class DiffusionOperator(TimeDependentOperator, LaplacianOperatorBase): def __init__(self, dimensions, diffusion_tensor=None, - dirichlet_bc=grudge.data.make_tdep_constant(0), dirichlet_tag="dirichlet", - neumann_bc=grudge.data.make_tdep_constant(0), neumann_tag="neumann", + dirichlet_bc=0, dirichlet_tag="dirichlet", + neumann_bc=0, neumann_tag="neumann", scheme=CentralSecondDerivative()): self.dimensions = dimensions @@ -61,11 +60,14 @@ class DiffusionOperator(TimeDependentOperator, LaplacianOperatorBase): assert self.dimensions == discr.dimensions - from grudge.mesh import check_bc_coverage + #from grudge.mesh import check_bc_coverage check_bc_coverage(discr.mesh, [self.dirichlet_tag, self.neumann_tag]) return BoundDiffusionOperator(self, discr) + def m_sym_operator(self): + return self.sym_operator(apply_minv=True) + def estimate_timestep(self, discr, stepper=None, stepper_class=None, stepper_args=None, t=None, fields=None): @@ -85,7 +87,8 @@ class DiffusionOperator(TimeDependentOperator, LaplacianOperatorBase): -class BoundDiffusionOperator(grudge.iterative.OperatorBase): +#class BoundDiffusionOperator(grudge.iterative.OperatorBase): +class BoundDiffusionOperator(): """Returned by :meth:`DiffusionOperator.bind`.""" def __init__(self, diffusion_op, discr): @@ -105,6 +108,8 @@ class BoundDiffusionOperator(grudge.iterative.OperatorBase): self.poincare_mean_value_hack = ( len(self.discr.get_boundary(BTAG_ALL).nodes) == len(self.discr.get_boundary(diffusion_op.neumann_tag).nodes)) + + def __call__(self, t, u): dop = self.diffusion_op diff --git a/grudge/models/gas_dynamics/__init__.py b/grudge/models/gas_dynamics/__init__.py index d50653c1..8afac735 100644 --- a/grudge/models/gas_dynamics/__init__.py +++ b/grudge/models/gas_dynamics/__init__.py @@ -763,7 +763,7 @@ class GasDynamicsOperator(TimeDependentOperator): #for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) return (fluxy(sym.int_tpair(int_operand, self.quad_tag)) + sum( - fluxy(sym.bv_tpair(sym.DOFDesc(sym.DTAG_BOUNDARY(tag), self.quad_tag),int_operand,join_fields(0, bdry_state, *bdry_fluxes) )) + fluxy(sym.bv_tpair(sym.DOFDesc(sym.DTAG_BOUNDARY(tag), self.quad_tag),int_operand,join_fields(1, bdry_state, *bdry_fluxes) )) for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) diff --git a/grudge/models/poisson.py b/grudge/models/poisson.py index a9f71457..81e4b56a 100644 --- a/grudge/models/poisson.py +++ b/grudge/models/poisson.py @@ -30,9 +30,9 @@ THE SOFTWARE. import numpy as np from grudge.models import Operator -from grudge.second_order import LDGSecondDerivative -import grudge.data -import grudge.iterative +from grudge.models.second_order import LDGSecondDerivative +from grudge import sym +#import grudge.iterative class LaplacianOperatorBase(object): @@ -48,15 +48,14 @@ class LaplacianOperatorBase(object): :class:`grudge.models.diffusion.DiffusionOperator` needs this. """ - from grudge.symbolic import Field, make_sym_vector - from grudge.second_order import SecondDerivativeTarget + from grudge.models.second_order import SecondDerivativeTarget if u is None: - u = Field("u") + u = sym.var("u") if dir_bc is None: - dir_bc = Field("dir_bc") + dir_bc = sym.var("dir_bc") if neu_bc is None: - neu_bc = Field("neu_bc") + neu_bc = sym.var("neu_bc") # strong_form here allows IPDG to reuse the value of grad u. grad_tgt = SecondDerivativeTarget( @@ -66,10 +65,11 @@ class LaplacianOperatorBase(object): def grad_bc_getter(tag, expr): assert tag == self.dirichlet_tag return dir_bc - self.scheme.grad(grad_tgt, + + self.scheme.grad(u, grad_tgt, bc_getter=grad_bc_getter, dirichlet_tags=[self.dirichlet_tag], - neumann_tags=[self.neumann_tag]) + neumann_tags=[self.neumann_tag] ) def apply_diff_tensor(v): if isinstance(self.diffusion_tensor, np.ndarray): @@ -115,9 +115,9 @@ class PoissonOperator(Operator, LaplacianOperatorBase): """ def __init__(self, dimensions, diffusion_tensor=None, - dirichlet_bc=grudge.data.ConstantGivenFunction(), + dirichlet_bc=0, dirichlet_tag="dirichlet", - neumann_bc=grudge.data.ConstantGivenFunction(), + neumann_bc=0, neumann_tag="neumann", scheme=LDGSecondDerivative()): self.dimensions = dimensions @@ -145,11 +145,12 @@ class PoissonOperator(Operator, LaplacianOperatorBase): return BoundPoissonOperator(self, discr) -class BoundPoissonOperator(grudge.iterative.OperatorBase): +#class BoundPoissonOperator(grudge.iterative.OperatorBase): +class BoundPoissonOperator(): """Returned by :meth:`PoissonOperator.bind`.""" def __init__(self, poisson_op, discr): - grudge.iterative.OperatorBase.__init__(self) + #grudge.iterative.OperatorBase.__init__(self) self.discr = discr pop = self.poisson_op = poisson_op diff --git a/grudge/models/second_order.py b/grudge/models/second_order.py index 24260c0f..b62cc988 100644 --- a/grudge/models/second_order.py +++ b/grudge/models/second_order.py @@ -257,9 +257,12 @@ class SecondDerivativeTarget(object): else: return vec*operand - def normal_times_flux(self, flux): - from grudge.flux import make_normal - return self.vec_times(make_normal(self.dimensions), flux) + def normal_times_flux(self, flux, dd=None): + + if dd is None: + return self.vec_times(sym.normal(flux.dd, self.dimensions), flux) + else: + return self.vec_times(sym.normal(dd, self.dimensions), flux) def apply_diff(self, nabla, operand): from pytools.obj_array import make_obj_array, is_obj_array @@ -275,11 +278,9 @@ class SecondDerivativeTarget(object): def _local_nabla(self): if self.strong_form: - from grudge.symbolic import make_stiffness - return make_stiffness(self.dimensions) + return sym.stiffness(self.dimensions) else: - from grudge.symbolic import make_stiffness_t - return make_stiffness_t(self.dimensions) + return sym.stiffness_t(self.dimensions) def add_derivative(self, operand=None): if operand is None: @@ -288,10 +289,15 @@ class SecondDerivativeTarget(object): self.add_local_derivatives( self.apply_diff(self._local_nabla(), operand)) - def add_inner_fluxes(self, flux, expr): - from grudge.symbolic import get_flux_operator + def flux_face_interp(self, pair): + return sym.interp(pair.dd, "all_faces")(pair) + + def add_inner_fluxes(self, flux): + #from grudge.symbolic import get_flux_operator + #self.inner_fluxes = self.inner_fluxes \ + #+ get_flux_operator(self.strong_neg*flux)(expr) self.inner_fluxes = self.inner_fluxes \ - + get_flux_operator(self.strong_neg*flux)(expr) + + self.flux_face_interp(self.strong_neg*flux) def add_boundary_flux(self, flux, volume_expr, bdry_expr, tag): from grudge.symbolic import BoundaryPair, get_flux_operator @@ -309,10 +315,9 @@ class SecondDerivativeTarget(object): @property def minv_all(self): - from grudge.symbolic.primitives import make_common_subexpression as cse from grudge.symbolic.operators import InverseMassOperator - return (cse(InverseMassOperator()(self.local_derivatives), "grad_loc") - + cse(InverseMassOperator()(self.fluxes), "grad_flux")) + return (sym.cse(InverseMassOperator()(self.local_derivatives), "grad_loc") + + sym.cse(InverseMassOperator()(self.fluxes), "grad_flux")) # }}} @@ -320,7 +325,7 @@ class SecondDerivativeTarget(object): # {{{ second derivative schemes class SecondDerivativeBase(object): - def grad(self, tgt, bc_getter, dirichlet_tags, neumann_tags): + def grad(self,u, tgt, bc_getter, dirichlet_tags, neumann_tags): """ :param bc_getter: a function (tag, volume_expr) -> boundary expr. *volume_expr* will be None to query the Neumann condition. @@ -328,19 +333,16 @@ class SecondDerivativeBase(object): n_times = tgt.normal_times_flux if tgt.strong_form: - def adjust_flux(f): - return n_times(u.int) - f + def adjust_flux(u, f,dd): + return n_times(u.int ,dd) - f else: - def adjust_flux(f): + def adjust_flux(u, f, dd): return f - - from grudge.flux import FluxScalarPlaceholder - u = FluxScalarPlaceholder() + int_u = sym.int_tpair(u) tgt.add_derivative() tgt.add_inner_fluxes( - adjust_flux(self.grad_interior_flux(tgt, u)), - tgt.int_flux_operand) + adjust_flux(int_u, self.grad_interior_flux(tgt, int_u), int_u.dd)) for tag in dirichlet_tags: tgt.add_boundary_flux( @@ -395,13 +397,12 @@ class LDGSecondDerivative(SecondDerivativeBase): return np.array([self.beta_value]*tgt.dimensions, dtype=np.float64) def grad_interior_flux(self, tgt, u): - from grudge.symbolic.primitives import make_common_subexpression as cse n_times = tgt.normal_times_flux v_times = tgt.vec_times return n_times( - cse(u.avg, "u_avg") - - v_times(self.beta(tgt), n_times(u.int-u.ext))) + sym.cse(u.avg, "u_avg") + - v_times(self.beta(tgt), n_times(u.int-u.ext, u.dd)), u.dd) def div(self, tgt, bc_getter, dirichlet_tags, neumann_tags): """ @@ -409,7 +410,6 @@ class LDGSecondDerivative(SecondDerivativeBase): *volume_expr* will be None to query the Neumann condition. """ - from grudge.symbolic.primitives import make_common_subexpression as cse from grudge.flux import FluxVectorPlaceholder, PenaltyTerm n_times = tgt.normal_times_flux @@ -431,13 +431,13 @@ class LDGSecondDerivative(SecondDerivativeBase): flux = n_times(flux_v.avg + v_times(self.beta(tgt), - cse(n_times(flux_v.int - flux_v.ext), "jump_v")) + sym.cse(n_times(flux_v.int - flux_v.ext), "jump_v")) - stab_term) from pytools.obj_array import make_obj_array - flux_arg_int = cse(make_obj_array(stab_term_generator.flux_args)) + flux_arg_int = sym.cse(make_obj_array(stab_term_generator.flux_args)) - tgt.add_derivative(cse(tgt.operand)) + tgt.add_derivative(sym.cse(tgt.operand)) tgt.add_inner_fluxes(adjust_flux(flux), flux_arg_int) self.add_div_bcs(tgt, bc_getter, dirichlet_tags, neumann_tags, @@ -459,9 +459,8 @@ class IPDGSecondDerivative(SecondDerivativeBase): self.stab_coefficient = stab_coefficient def grad_interior_flux(self, tgt, u): - from grudge.symbolic.primitives import make_common_subexpression as cse n_times = tgt.normal_times_flux - return n_times(cse(u.avg, "u_avg")) + return n_times(sym.cse(u.avg, "u_avg")) def div(self, tgt, bc_getter, dirichlet_tags, neumann_tags): """ @@ -469,7 +468,6 @@ class IPDGSecondDerivative(SecondDerivativeBase): *volume_expr* will be None to query the Neumann condition. """ - from grudge.symbolic.primitives import make_common_subexpression as cse from grudge.flux import FluxVectorPlaceholder, PenaltyTerm n_times = tgt.normal_times_flux @@ -496,9 +494,9 @@ class IPDGSecondDerivative(SecondDerivativeBase): flux = n_times(pure_diff_v.avg - stab_term) from pytools.obj_array import make_obj_array - flux_arg_int = cse(make_obj_array(stab_term_generator.flux_args)) + flux_arg_int = sym.cse(make_obj_array(stab_term_generator.flux_args)) - tgt.add_derivative(cse(tgt.operand)) + tgt.add_derivative(sym.cse(tgt.operand)) tgt.add_inner_fluxes(adjust_flux(flux), flux_arg_int) self.add_div_bcs(tgt, bc_getter, dirichlet_tags, neumann_tags, -- GitLab From 2a56429ddedd94793d5346a83ba6d292d25973e0 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Wed, 13 Jun 2018 00:25:47 -0500 Subject: [PATCH 12/15] When in doubt, destroy everything and start from advection --- examples/gas_dynamics/euler/vortex.py | 54 ++- examples/gas_dynamics/euler/vortex2.py | 175 +++++++++ .../gas_dynamics/gas_dynamics_initials.py | 52 +-- grudge/execution.py | 10 +- grudge/models/euler2.py | 178 +++++++++ grudge/models/gas_dynamics/__init__.py | 367 ++---------------- grudge/symbolic/dofdesc_inference.py | 2 + grudge/symbolic/operators.py | 4 +- grudge/symbolic/primitives.py | 4 +- 9 files changed, 431 insertions(+), 415 deletions(-) create mode 100644 examples/gas_dynamics/euler/vortex2.py create mode 100644 grudge/models/euler2.py diff --git a/examples/gas_dynamics/euler/vortex.py b/examples/gas_dynamics/euler/vortex.py index 7131f1b5..dd11a516 100644 --- a/examples/gas_dynamics/euler/vortex.py +++ b/examples/gas_dynamics/euler/vortex.py @@ -40,7 +40,7 @@ def main(write_output=True, visualize=False): for order in [3, 4, 5]: refine = 4 - mesh = generate_regular_rect_mesh(a=(0,-5), b=(10,5), n=(4,4), order=order) + mesh = generate_regular_rect_mesh(a=(0,-5), b=(10,5), n=(6,6), order=order) from gas_dynamics_initials import Vortex flow = Vortex() @@ -67,14 +67,18 @@ def main(write_output=True, visualize=False): bound_flow_vol = bind(discr, flow.volume_interpolant()) - fields = bound_flow_vol(queue, t=0) + fields = bound_flow_vol(queue, t=1) print(fields) euler_ex = bind(discr, op.sym_operator()) max_eigval = [0] def rhs(t, q): - ode_rhs, speed = euler_ex(t, q) + #ode_rhs, speed = euler_ex(queue, t=t, q=q) + ting = euler_ex(queue, t=t, q=q) + ode_rhs = ting[:-1] + speed = ting[-1] + #pu.db max_eigval[0] = speed return ode_rhs rhs(0, fields) @@ -89,19 +93,34 @@ def main(write_output=True, visualize=False): #stepper = SSP3TimeStepper( #vector_primitive_factory=discr.get_vector_primitive_factory()) from grudge.shortcuts import set_up_rk4 - dt = 0.000001 - final_t = 0.0001 + dt = 0.001 + final_t = 0.1 dt_stepper = set_up_rk4("q", dt, fields, rhs) if visualize: from grudge.shortcuts import make_visualizer vis = make_visualizer(discr, vis_order=order) - step = 0 + norm = bind(discr, sym.norm(2, sym.var("q"))) + + vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, 0), + [ + ("rho", op.rho(fields)), + ("e", op.e(fields)), + ("rho_u", op.rho_u(fields)), + ("u", op.u(fields)) + #("true_rho", op.rho(true_fields)), + #("true_e", op.e(true_fields)), + #("true_rho_u", op.rho_u(true_fields)), + #("true_u", op.u(true_fields) + #("rhs_rho", op.rho(rhs_fields)), + #("rhs_e", op.e(rhs_fields)), + #("rhs_rho_u", op.rho_u(rhs_fields)), + ]) + step = 1 for event in dt_stepper.run(t_end=final_t): if isinstance(event, dt_stepper.StateComputed): - #assert event.component_id == "u" - print(event.component_id) + assert event.component_id == "q" fields = event.state_component if visualize: @@ -125,17 +144,16 @@ def main(write_output=True, visualize=False): print(step) - pu.db - true_fields = flow.volume_interpolant(final_time, discr) - l2_error = discr.norm(fields-true_fields) - l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields)) - l2_error_e = discr.norm(op.e(fields)-op.e(true_fields)) - l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields)) - l2_error_u = discr.norm(op.u(fields)-op.u(true_fields)) + #true_fields = bound_flow_vol(queue, t=final_t) + #l2_error = norm(queue, q = fields-true_fields) + #l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields)) + #l2_error_e = discr.norm(op.e(fields)-op.e(true_fields)) + #l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields)) + #l2_error_u = discr.norm(op.u(fields)-op.u(true_fields)) - eoc_rec.add_data_point(order, l2_error) - print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) + #eoc_rec.add_data_point(order, l2_error) + #print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) @@ -146,7 +164,7 @@ def main(write_output=True, visualize=False): if __name__ == "__main__": - main() + main(visualize=True) diff --git a/examples/gas_dynamics/euler/vortex2.py b/examples/gas_dynamics/euler/vortex2.py new file mode 100644 index 00000000..4a4871cd --- /dev/null +++ b/examples/gas_dynamics/euler/vortex2.py @@ -0,0 +1,175 @@ +# Copyright (C) 2008 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 . + + + + +from __future__ import division +from __future__ import absolute_import +from __future__ import print_function +import numpy + +import pyopencl as cl # noqa + +from grudge import sym, bind, DGDiscretizationWithBoundaries + +def main(write_output=True, visualize=False): + from pytools import add_python_path_relative_to_script + add_python_path_relative_to_script("..") + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + from meshmode.mesh.generation import generate_regular_rect_mesh + + + for order in [3, 4, 5]: + mesh = generate_regular_rect_mesh(a=(0,-5), b=(10,5), n=(6,6), order=order) + from gas_dynamics_initials import Vortex + flow = Vortex() + + from grudge.models.euler2 import EulerOperator + from grudge.models.gas_dynamics import ( + GasDynamicsOperator, PolytropeEOS, GammaLawEOS) + + from meshmode.mesh import BTAG_ALL + from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory + # works equally well for GammaLawEOS + op = EulerOperator(dimensions=2, mu=flow.mu, + + equation_of_state=PolytropeEOS(flow.gamma), + bc_inflow=flow.volume_interpolant(), + inflow_tag=BTAG_ALL) + + discr = DGDiscretizationWithBoundaries(cl_ctx,mesh, order=order, + quad_tag_to_group_factory={ + "product": QuadratureSimplexGroupFactory(order=3*order) + }) + + if visualize: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + + + bound_flow_vol = bind(discr, flow.volume_interpolant()) + fields = bound_flow_vol(queue, t=1) + print(fields) + + euler_ex = bind(discr, op.sym_operator()) + + max_eigval = [0] + def rhs(t, q): + #ode_rhs, speed = euler_ex(queue, t=t, q=q) + ting = euler_ex(queue, t=t, q=q) + ode_rhs = ting[:-1] + speed = ting[-1] + #pu.db + max_eigval[0] = speed + return ode_rhs + rhs(0, fields) + + + # limiter ------------------------------------------------------------ + #from grudge.models.gas_dynamics import SlopeLimiter1NEuler + #limiter = SlopeLimiter1NEuler(discr, flow.gamma, 2, op) + + #from grudge.timestep.runge_kutta import SSP3TimeStepper + #stepper = SSP3TimeStepper(limiter=limiter) + #stepper = SSP3TimeStepper( + #vector_primitive_factory=discr.get_vector_primitive_factory()) + from grudge.shortcuts import set_up_rk4 + dt = 0.001 + final_t = 0.1 + dt_stepper = set_up_rk4("q", dt, fields, rhs) + + if visualize: + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=order) + + norm = bind(discr, sym.norm(2, sym.var("q"))) + + vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, 0), + [ + ("rho", op.rho(fields)), + ("e", op.e(fields)), + ("rho_u", op.rho_u(fields)), + ("u", op.u(fields)) + #("true_rho", op.rho(true_fields)), + #("true_e", op.e(true_fields)), + #("true_rho_u", op.rho_u(true_fields)), + #("true_u", op.u(true_fields) + #("rhs_rho", op.rho(rhs_fields)), + #("rhs_e", op.e(rhs_fields)), + #("rhs_rho_u", op.rho_u(rhs_fields)), + ]) + step = 1 + for event in dt_stepper.run(t_end=final_t): + if isinstance(event, dt_stepper.StateComputed): + assert event.component_id == "q" + fields = event.state_component + + if visualize: + vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, step), + [ + ("rho", op.rho(fields)), + ("e", op.e(fields)), + ("rho_u", op.rho_u(fields)), + ("u", op.u(fields)), + + #("true_rho", op.rho(true_fields)), + #("true_e", op.e(true_fields)), + #("true_rho_u", op.rho_u(true_fields)), + #("true_u", op.u(true_fields)), + + #("rhs_rho", op.rho(rhs_fields)), + #("rhs_e", op.e(rhs_fields)), + #("rhs_rho_u", op.rho_u(rhs_fields)), + ]) + step += 1 + print(step) + + + + #true_fields = bound_flow_vol(queue, t=final_t) + #l2_error = norm(queue, q = fields-true_fields) + #l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields)) + #l2_error_e = discr.norm(op.e(fields)-op.e(true_fields)) + #l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields)) + #l2_error_u = discr.norm(op.u(fields)-op.u(true_fields)) + + #eoc_rec.add_data_point(order, l2_error) + #print(eoc_rec.pretty_print("P.Deg.", "L2 Error")) + + + + # after order loop + assert eoc_rec.estimate_order_of_convergence()[0,1] > 6 + + + + +if __name__ == "__main__": + main(visualize=True) + + + +# entry points for py.test ---------------------------------------------------- +from pytools.test import mark_test +@mark_test.long +def test_euler_vortex(): + main(write_output=False) diff --git a/examples/gas_dynamics/gas_dynamics_initials.py b/examples/gas_dynamics/gas_dynamics_initials.py index de5db3d3..9d24aa7f 100644 --- a/examples/gas_dynamics/gas_dynamics_initials.py +++ b/examples/gas_dynamics/gas_dynamics_initials.py @@ -121,50 +121,6 @@ class UniformMachFlow: tag=tag, kind=discr.compute_kind, dtype=discr.default_scalar_type) -class Vortex: - def __init__(self): - self.beta = 5 - self.gamma = 1.4 - self.center = numpy.array([5, 0]) - self.velocity = numpy.array([1, 0]) - - self.mu = 0 - self.prandtl = 0.72 - self.spec_gas_const = 287.1 - - def __call__(self): - x_vec = sym.nodes(2) - t = sym.var("t", sym.DD_SCALAR) - - vortex_loc = self.center + t*self.velocity - - # coordinates relative to vortex center - x_rel = x_vec[0] - vortex_loc[0] - y_rel = x_vec[1] - vortex_loc[1] - - # Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159 - # also JSH/TW Nodal DG Methods, p. 209 - - from math import pi - r = numpy.sqrt(x_rel**2+y_rel**2) - expterm = self.beta*numpy.exp(1-r**2) - u = self.velocity[0] - expterm*y_rel/(2*pi) - v = self.velocity[1] + expterm*x_rel/(2*pi) - rho = (1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1)) - p = rho**self.gamma - - e = p/(self.gamma-1) + rho/2*(u**2+v**2) - - from grudge.tools import join_fields - return join_fields(rho, e, rho*u, rho*v) - - - def boundary_interpolant(self, t, discr, tag): - return discr.convert_boundary( - self(t, discr.get_boundary(tag).nodes.T - .astype(discr.default_scalar_type)), - tag=tag, kind=discr.compute_kind) - @@ -206,17 +162,11 @@ class Vortex: e = p/(self.gamma-1) + rho/2*(u**2+v**2) from pytools.obj_array import join_fields - return join_fields(rho+100, e+100, rho*u+100, rho*v+100) + return join_fields(rho, e, rho*u, rho*v) def volume_interpolant(self): return self() - def boundary_interpolant(self, t, discr, tag): - return discr.convert_boundary( - self(t, discr.get_boundary(tag).nodes.T - .astype(discr.default_scalar_type)), - tag=tag, kind=discr.compute_kind) - diff --git a/grudge/execution.py b/grudge/execution.py index 6e942f43..b7e24d48 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -390,12 +390,12 @@ class ExecutionMapper(mappers.Evaluator, % (name, type(v))) kwargs["grdg_n"] = discr.nnodes - kdescr.loopy_kernel = lp.set_options(kdescr.loopy_kernel, trace_assignment_values=True) + #kdescr.loopy_kernel = lp.set_options(kdescr.loopy_kernel, trace_assignment_values=True) evt, result_dict = kdescr.loopy_kernel(self.queue, **kwargs) - if "expr_217" in result_dict.keys(): - knl = lp.set_options(kdescr.loopy_kernel, "write_cl") - knl(self.queue, **kwargs) - pu.db + #if "expr_217" in result_dict.keys(): + #knl = lp.set_options(kdescr.loopy_kernel, "write_cl") + #knl(self.queue, **kwargs) + #pu.db return list(result_dict.items()), [] def map_insn_assign(self, insn): diff --git a/grudge/models/euler2.py b/grudge/models/euler2.py new file mode 100644 index 00000000..a842497f --- /dev/null +++ b/grudge/models/euler2.py @@ -0,0 +1,178 @@ +# -*- coding: utf8 -*- +"""Operators modeling advective phenomena.""" + +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2009-2017 Andreas Kloeckner, Bogdan Enache" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + +import numpy as np +import numpy +import numpy.linalg as la + +from grudge.models import HyperbolicOperator +from grudge import sym +from pytools.obj_array import make_obj_array, join_fields + + + +class EquationOfState(object): + def q_to_p(self, op, q): + raise NotImplementedError + + def p_to_e(self, p, rho, u): + raise NotImplementedError + +class GammaLawEOS(EquationOfState): + # FIXME Shouldn't gamma only occur in the equation of state? + # I.e. shouldn't all uses of gamma go through the EOS? + + def __init__(self, gamma): + self.gamma = gamma + + def q_to_p(self, op, q): + return (self.gamma-1)*(op.e(q)-0.5*numpy.dot(op.rho_u(q), op.u(q))) + + def p_to_e(self, p, rho, u): + return p / (self.gamma - 1) + rho / 2 * numpy.dot(u, u) + +class PolytropeEOS(GammaLawEOS): + # inverse is same as superclass + + def q_to_p(self, op, q): + return op.rho(q)**self.gamma + + + + + +class EulerOperator(HyperbolicOperator): + def __init__(self, dimensions, + gamma=None, mu=0, bc_inflow=None, + bc_outflow=None, + bc_noslip=None, + bc_supersonic_inflow=None, + prandtl=None, spec_gas_const=1.0, + equation_of_state=None, + inflow_tag="inflow", + outflow_tag="outflow", + noslip_tag="noslip", + wall_tag="wall", + supersonic_inflow_tag="supersonic_inflow", + supersonic_outflow_tag="supersonic_outflow", + source=None, + second_order_scheme=None, + artificial_viscosity_mode=None, + quad_tag = "product" + ): + self.dimensions = dimensions + self.equation_of_state = equation_of_state + + + def rho(self, q): + return q[0] + + def e(self, q): + return q[1] + + def rho_u(self, q): + return q[2:2+self.dimensions] + + def u(self, q): + rho_u = self.rho_u(q) + return make_obj_array([ + rho_u[i]/self.rho(q) + for i in range(self.dimensions)]) + + def p(self,q): + return self.equation_of_state.q_to_p(self, q) + + def flux(self, q): + from pytools import delta + + return [ # one entry for each flux direction + join_fields( + # flux rho + self.rho_u(q)[i], + + # flux E + (self.e(q)+self.p(q))*self.u(q)[i], + + # flux rho_u + make_obj_array([ + self.rho_u(q)[i]*self.u(q)[j] + + delta(i,j) * self.p(q) + for j in range(self.dimensions) + ]) + ) + for i in range(self.dimensions)] + + + def lax_flux(self, speed, state, flux): + + normal = sym.normal(state.dd, self.dimensions) + + wave_speed_ph = speed + state_ph = state + fluxes_ph = flux + + from pymbolic.primitives import Max + #penalty = Max((wave_speed_ph.int,wave_speed_ph.ext))*(state_ph.ext-state_ph.int) + penalty = 0 + + num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) + - penalty) + return num_flux + + def characteristic_velocity(self,state): + from grudge.symbolic.operators import ElementwiseMaxOperator + from grudge.symbolic.primitives import CFunction + sqrt = CFunction("sqrt") + + sound_speed = sqrt(self.equation_of_state.gamma*self.p(state)/self.rho(state)) + u = self.u(state) + speed = sqrt(numpy.dot(u, u)) + sound_speed + return ElementwiseMaxOperator("vol","vol")(speed) + + + + def sym_operator(self): + q = sym.make_sym_array("q", self.dimensions+2) + + + def to_all_faces(pair, dd): + return sym.interp(dd, "all_faces")(self.flux(pair)) + + speed = 0 #self.characteristic_velocity(q) + volq_flux = self.flux(q) + stiff_t = sym.stiffness_t(self.dimensions) + ans = sym.InverseMassOperator()( + np.dot(stiff_t, volq_flux) + - sym.FaceMassOperator()( + to_all_faces(self.lax_flux(sym.int_tpair(speed), sym.int_tpair(q), sym.int_tpair(self.flux(q)) ), sym.int_tpair(speed).dd) + #+ flux(sym.bv_tpair(sym.BTAG_ALL, q, bc_in)) + + ) + ) + return join_fields(ans,speed) +# vim: foldmethod=marker diff --git a/grudge/models/gas_dynamics/__init__.py b/grudge/models/gas_dynamics/__init__.py index 8afac735..34783833 100644 --- a/grudge/models/gas_dynamics/__init__.py +++ b/grudge/models/gas_dynamics/__init__.py @@ -36,16 +36,15 @@ import meshmode.mesh from grudge.models import TimeDependentOperator from pytools import Record from grudge.tools import is_zero -from pymbolic.primitives import make_common_subexpression as cse +#from pymbolic.primitives import make_common_subexpression as cse from pytools import memoize_method from pytools.obj_array import make_obj_array, join_fields from grudge import sym AXES = ["x", "y", "z", "w"] - - - +def cse(x, stsr=None): + return x # {{{ equations of state @@ -78,9 +77,6 @@ class PolytropeEOS(GammaLawEOS): # }}} - - - class GasDynamicsOperator(TimeDependentOperator): """An nD Navier-Stokes and Euler operator. @@ -147,7 +143,10 @@ class GasDynamicsOperator(TimeDependentOperator): equation_of_state = GammaLawEOS(gamma) - dull_bc = make_obj_array( [1, 1] + [0]*dimensions) + mostly_one = sym.nodes(dimensions)/sym.nodes(dimensions) + mostly_zero = sym.nodes(dimensions)*0 + #dull_bc = make_obj_array( [1, 1] + [0]*dimensions) + dull_bc = make_obj_array( [mostly_one,mostly_one] + [mostly_one]*dimensions) if bc_inflow is None: bc_inflow = dull_bc if bc_outflow is None: @@ -187,7 +186,7 @@ class GasDynamicsOperator(TimeDependentOperator): raise ValueError("artificial_viscosity_mode has an invalid value") self.artificial_viscosity_mode = artificial_viscosity_mode - self.quad_tag = quad_tag + self.quad_tag = None #quad_tag @@ -200,6 +199,7 @@ class GasDynamicsOperator(TimeDependentOperator): @memoize_method def volq_state(self): + return self.state() quad_dd = sym.DOFDesc("vol", self.quad_tag) to_vol_quad = sym.interp("vol", quad_dd) return cse(to_vol_quad(self.state()), "vol_quad_state") @@ -245,12 +245,6 @@ class GasDynamicsOperator(TimeDependentOperator): def cse_p(self, q): return cse(self.p(q), "p") - def temperature(self, q): - c_v = 1 / (self.equation_of_state.gamma - 1) * self.spec_gas_const - return (self.e(q)/self.rho(q) - 0.5 * numpy.dot(self.u(q), self.u(q))) / c_v - - def cse_temperature(self, q): - return cse(self.temperature(q), "temperature") def get_mu(self, q, to_quad_op): """ @@ -318,10 +312,10 @@ class GasDynamicsOperator(TimeDependentOperator): from grudge.symbolic.primitives import CFunction sqrt = CFunction("sqrt") - #sound_speed = cse(sqrt( - #self.equation_of_state.gamma*self.cse_p(state)/self.cse_rho(state)), - #"sound_speed") - sound_speed = sqrt(self.cse_p(state)) + sound_speed = cse(sqrt( + self.equation_of_state.gamma*self.cse_p(state)/self.cse_rho(state)), + "sound_speed") + #sound_speed = sqrt(self.cse_p(state)) #sound_speed = self.cse_p(state) u = self.cse_u(state) speed = cse(sqrt(numpy.dot(u, u)), "norm_u") + sound_speed @@ -338,116 +332,7 @@ class GasDynamicsOperator(TimeDependentOperator): return do - # }}} - - # {{{ helpers for second-order part --------------------------------------- - - # {{{ compute gradient of state --------------------------------------- - def grad_of(self, var, faceq_var): - from grudge.second_order import SecondDerivativeTarget - grad_tgt = SecondDerivativeTarget( - self.dimensions, strong_form=False, - operand=var, - int_flux_operand=faceq_var, - bdry_flux_int_operand=faceq_var) - - #self.second_order_scheme.grad(grad_tgt, - #bc_getter=self.get_boundary_condition_for, - #dirichlet_tags=self.get_boundary_tags(), - #neumann_tags=[]) - - return grad_tgt.minv_all - - def grad_of_state(self): - dimensions = self.dimensions - - state = self.state() - - dq = numpy.zeros((len(state), dimensions), dtype=object) - - for i in range(len(state)): - dq[i,:] = self.grad_of( - state[i], self.faceq_state()[i]) - - return dq - - def grad_of_state_func(self, func, of_what_descr): - return cse(self.grad_of( - func(self.volq_state()), - func(self.faceq_state())), - "grad_"+of_what_descr) - - # }}} - - # {{{ viscous stress tensor - - def tau(self, to_quad_op, state, mu=None): - #faceq_state = self.faceq_state() - - dimensions = self.dimensions - - # {{{ compute gradient of u --------------------------------------- - # Use the product rule to compute the gradient of - # u from the gradient of (rho u). This ensures we don't - # compute the derivatives twice. - - from pytools.obj_array import with_object_array_or_scalar - dq = with_object_array_or_scalar( - to_quad_op, self.grad_of_state()) - - q = cse(to_quad_op(state)) - - du = numpy.zeros((dimensions, dimensions), dtype=object) - for i in range(dimensions): - for j in range(dimensions): - du[i,j] = cse( - (dq[i+2,j] - self.cse_u(q)[i] * dq[0,j]) / self.rho(q), - "du%d_d%s" % (i, AXES[j])) - - # }}} - - # {{{ put together viscous stress tau ----------------------------- - from pytools import delta - - if mu is None: - mu = self.get_mu(q, to_quad_op) - - tau = numpy.zeros((dimensions, dimensions), dtype=object) - for i in range(dimensions): - for j in range(dimensions): - tau[i,j] = cse(mu * cse(du[i,j] + du[j,i] - - 2/self.dimensions * delta(i,j) * numpy.trace(du)), - "tau_%d%d" % (i, j)) - - return tau - - # }}} - - # }}} - - # }}} - - # {{{ heat conduction - - def heat_conduction_coefficient(self, to_quad_op): - mu = self.get_mu(self.state(), to_quad_op) - if self.prandtl is None or numpy.isinf(self.prandtl): - return 0 - - eos = self.equation_of_state - return (mu / self.prandtl) * (eos.gamma / (eos.gamma-1)) - - def heat_conduction_grad(self, to_quad_op): - grad_p_over_rho = self.grad_of_state_func( - lambda state: self.p(state)/self.rho(state), - "p_over_rho") - return (self.heat_conduction_coefficient(to_quad_op) - * to_quad_op(grad_p_over_rho)) - - # }}} - - # {{{ flux def flux(self, q): from pytools import delta @@ -479,9 +364,10 @@ class GasDynamicsOperator(TimeDependentOperator): BC is linearized. """ if state0 is None: - state0 = sym.make_sym_array(bc_name, self.dimensions+2) + raise ValueError + #state0 = sym.make_sym_array(bc_name, self.dimensions+2) - boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) + #boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) boundary_dd = sym.DOFDesc(sym.DTAG_BOUNDARY(tag), self.quad_tag) to_bdry_quad = sym.interp(sym.BTAG_ALL, boundary_dd) #state0 = cse(to_bdry_quad(state0)) @@ -514,7 +400,8 @@ class GasDynamicsOperator(TimeDependentOperator): def outflow_state(self, state): normal = sym.normal(sym.DTAG_BOUNDARY(self.outflow_tag), self.dimensions) - bc = self.make_bc_info("bc_q_out", self.outflow_tag, state) + #bc = self.make_bc_info("bc_q_out", self.outflow_tag, state, self.bc_outflow) + bc = self.make_bc_info("bc_q_out", self.outflow_tag, state, state) # see grudge/doc/maxima/euler.mac return join_fields( @@ -549,15 +436,14 @@ class GasDynamicsOperator(TimeDependentOperator): def inflow_state(self, state): normal = sym.normal(sym.DTAG_BOUNDARY(self.inflow_tag), self.dimensions) - bc = self.make_bc_info("bc_q_in", self.inflow_tag, state) + #bc = self.make_bc_info("bc_q_in", self.inflow_tag, state, self.bc_inflow) + bc = self.make_bc_info("bc_q_in", self.inflow_tag, state, state) return self.inflow_state_inner(normal, bc, "inflow") def noslip_state(self, state): - state0 = join_fields( - sym.make_sym_array("bc_q_noslip", 2), - [0]*self.dimensions) normal = sym.normal(sym.DTAG_BOUNDARY(self.noslip_tag), self.dimensions) - bc = self.make_bc_info("bc_q_noslip", self.noslip_tag, state, state0) + #bc = self.make_bc_info("bc_q_noslip", self.noslip_tag, state, self.bc_noslip) + bc = self.make_bc_info("bc_q_noslip", self.noslip_tag, state, state) return self.inflow_state_inner(normal, bc, "noslip") def wall_state(self, state): @@ -589,8 +475,8 @@ class GasDynamicsOperator(TimeDependentOperator): state = self.state() return { - self.supersonic_inflow_tag: - sym.make_sym_array("bc_q_supersonic_in", self.dimensions+2), + self.supersonic_inflow_tag: self.bc_supersonic_inflow, + #sym.make_sym_array("bc_q_supersonic_in", self.dimensions+2), self.supersonic_outflow_tag: #sym.interp("vol",self.supersonic_outflow_tag)(state), #TODO: It should be this according to my em translation <- but (it spits out an error) sym.interp("vol",sym.DTAG_BOUNDARY(self.supersonic_outflow_tag))(state), #TODO: Attempt 2 @@ -600,6 +486,7 @@ class GasDynamicsOperator(TimeDependentOperator): @memoize_method def get_boundary_tags(self): + return set([self.inflow_tag,]) return (set(self.get_primitive_boundary_conditions().keys()) | set(self.get_conservative_boundary_conditions().keys())) @@ -627,7 +514,7 @@ class GasDynamicsOperator(TimeDependentOperator): to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd) if tag in prim_bcs: - # BC is given in primitive variables, avoid converting + # BC is giVEn in primitive variables, avoid converting # to conservative and back. try: norm_expr = self._normalize_expr(expr) @@ -665,51 +552,9 @@ class GasDynamicsOperator(TimeDependentOperator): # }}} - # {{{ second order part - def div(self, vol_operand, int_face_operand): - from grudge.second_order import SecondDerivativeTarget - div_tgt = SecondDerivativeTarget( - self.dimensions, strong_form=False, - operand=vol_operand, - int_flux_operand=int_face_operand) - - #self.second_order_scheme.div(div_tgt, - #bc_getter=self.get_boundary_condition_for, - #dirichlet_tags=list(self.get_boundary_tags()), - #neumann_tags=[]) - - return div_tgt.minv_all - - #def make_second_order_part(self): - #state = self.state() - #faceq_state = self.faceq_state() - #volq_state = self.volq_state() - #int_faces_dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR, self.quad_tag) - #to_int_face_quad = sym.interp("vol", int_faces_dd) -# - #volq_tau_mat = self.tau(to_vol_quad, state) - #faceq_tau_mat = self.tau(to_int_face_quad, state) -# - #return join_fields( - #0, - #self.div( - #numpy.sum(volq_tau_mat*self.cse_u(volq_state), axis=1) - #+ self.heat_conduction_grad(to_vol_quad) - #, - #numpy.sum(faceq_tau_mat*self.cse_u(faceq_state), axis=1) - #+ self.heat_conduction_grad(to_int_face_quad) - #, - #), - #[ - #self.div(volq_tau_mat[i], faceq_tau_mat[i]) - #for i in range(self.dimensions)] - #) - # }}} # {{{ operator template --------------------------------------------------- - def make_extra_terms(self): - return 0 def make_lax_friedrichs_flux(self, wave_speed, state, fluxes, bdry_tags_states_and_fluxes,strong): #from hedge.flux import FluxVectorPlaceholder, flux_max @@ -726,7 +571,8 @@ class GasDynamicsOperator(TimeDependentOperator): print(fvph.dd) from pymbolic.primitives import Max - penalty = Max((wave_speed_ph.int,wave_speed_ph.ext))*(state_ph.ext-state_ph.int) + #penalty = Max((wave_speed_ph.int,wave_speed_ph.ext))*(state_ph.ext-state_ph.int) + penalty = 0 if not strong: num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) # TODO dd error between flux and norm here @@ -740,28 +586,10 @@ class GasDynamicsOperator(TimeDependentOperator): def fluxy(pair): return sym.interp(pair.dd, all_faces_dd)(lff(pair)) - #def no_int_int_pair(expression): - #from grudge.symbolic.primitives import TracePair - #i = expression - #e = cse(sym.OppositeInteriorFaceSwap()(i)) - #return TracePair(sym.DOFDesc("int_faces", self.quad_tag), i, e) - - - - #from hedge.optemplate import get_flux_operator - #flux_op = get_flux_operator(num_flux) int_operand = join_fields(wave_speed, state, *fluxes) - - #return (flux(int_operand) + sum( - #flux(sym.bv_tpair(sym.DTAG_BOUNDARY(tag),int_operand,join_fields(0, bdry_state, *bdry_fluxes) )) - #for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) - - #return (flux(sym.int_tpair(int_operand, self.quad_tag)) + sum( - #flux(sym.bv_tpair(sym.DTAG_BOUNDARY(tag),int_operand,join_fields(0, bdry_state, *bdry_fluxes) )) - #for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) - + return fluxy(sym.int_tpair(int_operand, self.quad_tag)) return (fluxy(sym.int_tpair(int_operand, self.quad_tag)) + sum( fluxy(sym.bv_tpair(sym.DOFDesc(sym.DTAG_BOUNDARY(tag), self.quad_tag),int_operand,join_fields(1, bdry_state, *bdry_fluxes) )) for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) @@ -777,21 +605,6 @@ class GasDynamicsOperator(TimeDependentOperator): int_faces_dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR, self.quad_tag) to_int_face_quad = sym.interp("vol", int_faces_dd) - # {{{ artificial diffusion - #def make_artificial_diffusion(): - #if self.artificial_viscosity_mode not in ["diffusion"]: - #return 0 - - #dq = self.grad_of_state() - - #return make_obj_array([ - #self.div( - #to_vol_quad(self.sensor())*to_vol_quad(dq[i]), - #to_int_face_quad(self.sensor())*to_int_face_quad(dq[i])) - #for i in range(dq.shape[0])]) - # }}} - - # {{{ state setup volq_flux = self.flux(self.volq_state()) faceq_flux = self.flux(self.faceq_state()) @@ -852,125 +665,5 @@ class GasDynamicsOperator(TimeDependentOperator): return result - # }}} - - # }}} - - # {{{ operator binding ---------------------------------------------------- - def bind(self, discr, sensor=None, sensor_scaling=None, viscosity_only=False): - if (sensor is None and - self.artificial_viscosity_mode is not None): - raise ValueError("must specify a sensor if using " - "artificial viscosity") - - bound_op = discr.compile(self.sym_operator( - sensor_scaling=sensor_scaling, - viscosity_only=False)) - - from meshmode.mesh import check_bc_coverage - check_bc_coverage(discr.mesh, [ - self.inflow_tag, - self.outflow_tag, - self.noslip_tag, - self.wall_tag, - self.supersonic_inflow_tag, - self.supersonic_outflow_tag, - ]) - - if self.mu == 0 and not discr.get_boundary(self.noslip_tag).is_empty(): - raise RuntimeError("no-slip BCs only make sense for " - "viscous problems") - - def rhs(t, q): - extra_kwargs = {} - if self.source is not None: - extra_kwargs["source_vect"] = self.source.volume_interpolant( - t, q, discr) - - if sensor is not None: - extra_kwargs["sensor"] = sensor(q) - - opt_result = bound_op(q=q, - bc_q_in=self.bc_inflow.boundary_interpolant( - t, discr, self.inflow_tag), - bc_q_out=self.bc_outflow.boundary_interpolant( - t, discr, self.outflow_tag), - bc_q_noslip=self.bc_noslip.boundary_interpolant( - t, discr, self.noslip_tag), - bc_q_supersonic_in=self.bc_supersonic_inflow - .boundary_interpolant(t, discr, - self.supersonic_inflow_tag), - **extra_kwargs - ) - - max_speed = opt_result[-1] - ode_rhs = opt_result[:-1] - return ode_rhs, discr.nodewise_max(max_speed) - - return rhs - - # }}} - - # {{{ timestep estimation ------------------------------------------------- - - def estimate_timestep(self, discr, - stepper=None, stepper_class=None, stepper_args=None, - t=None, max_eigenvalue=None): - u"""Estimate the largest stable timestep, given a time stepper - `stepper_class`. If none is given, RK4 is assumed. - """ - - dg_factor = (discr.dt_non_geometric_factor() - * discr.dt_geometric_factor()) - - # see JSH/TW, eq. (7.32) - rk4_dt = dg_factor / (max_eigenvalue + self.mu / dg_factor) - - from grudge.timestep.stability import \ - approximate_rk4_relative_imag_stability_region - return rk4_dt * approximate_rk4_relative_imag_stability_region( - stepper, stepper_class, stepper_args) - - # }}} - - - - -# {{{ limiter (unfinished, deprecated) -class SlopeLimiter1NEuler: - def __init__(self, discr, gamma, dimensions, op): - """Construct a limiter from Jan's book page 225 - """ - self.discr = discr - self.gamma=gamma - self.dimensions=dimensions - self.op=op - - from grudge.symbolic.operators import AveragingOperator - self.get_average = AveragingOperator().bind(discr) - - def __call__(self, fields): - from grudge.tools import join_fields - - #get conserved fields - rho=self.op.rho(fields) - e=self.op.e(fields) - rho_velocity=self.op.rho_u(fields) - - #get primitive fields - #to do - - #reset field values to cell average - rhoLim=self.get_average(rho) - eLim=self.get_average(e) - temp=join_fields([self.get_average(rho_vel) - for rho_vel in rho_velocity]) - - #should do for primitive fields too - - return join_fields(rhoLim, eLim, temp) - -# }}} - - + # vim: foldmethod=marker diff --git a/grudge/symbolic/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py index 379969bc..942b805b 100644 --- a/grudge/symbolic/dofdesc_inference.py +++ b/grudge/symbolic/dofdesc_inference.py @@ -99,6 +99,8 @@ class DOFDescInferenceMapper(RecursiveMapper, CSECachingMapperMixin): self.name_to_dofdesc = name_to_dofdesc def infer_for_name(self, name): + #if(name=='expr_29'): + #pu.db try: return self.name_to_dofdesc[name] except KeyError: diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index fcf87d23..31dc3166 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -59,8 +59,8 @@ class Operator(pymbolic.primitives.Expression): self.dd_out = _sym().as_dofdesc(dd_out) #TODO: remove this. According to Python docs this creates weird circular references # And a bad time for everyone - import inspect - self.stack = inspect.stack() + #import inspect + #self.stack = inspect.stack() def stringifier(self): from grudge.symbolic.mappers import StringifyMapper diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 9fd7b790..a29e481a 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -201,8 +201,8 @@ class DOFDesc(object): #TODO: remove this. According to Python docs this creates weird circular references # And a bad time for everyone - import inspect - self.stack = inspect.stack() + #import inspect + #self.stack = inspect.stack() def is_scalar(self): -- GitLab From c61a3142d11ce3fdb56d14f942350e184a8b8fcf Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Thu, 14 Jun 2018 02:01:13 -0500 Subject: [PATCH 13/15] fixed a fluxy bug --- grudge/models/euler2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grudge/models/euler2.py b/grudge/models/euler2.py index a842497f..c7fb2224 100644 --- a/grudge/models/euler2.py +++ b/grudge/models/euler2.py @@ -161,7 +161,7 @@ class EulerOperator(HyperbolicOperator): def to_all_faces(pair, dd): - return sym.interp(dd, "all_faces")(self.flux(pair)) + return sym.interp(dd, "all_faces")(pair) speed = 0 #self.characteristic_velocity(q) volq_flux = self.flux(q) -- GitLab From 286b1e7cfcbcc1841e1e6456434d5580fa99d8d7 Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Tue, 19 Jun 2018 02:37:02 -0500 Subject: [PATCH 14/15] can advect things --- examples/advection/gas_dynamics_initials.py | 172 ++++++++++++++++++ examples/advection/new.py | 138 ++++++++++++++ examples/advection/var-velocity.py | 8 +- examples/gas_dynamics/euler/vortex2.py | 12 +- .../gas_dynamics/gas_dynamics_initials.py | 7 +- grudge/execution.py | 4 +- grudge/models/advection.py | 28 ++- grudge/models/euler2.py | 43 ++++- 8 files changed, 377 insertions(+), 35 deletions(-) create mode 100644 examples/advection/gas_dynamics_initials.py create mode 100644 examples/advection/new.py diff --git a/examples/advection/gas_dynamics_initials.py b/examples/advection/gas_dynamics_initials.py new file mode 100644 index 00000000..6f4e4f41 --- /dev/null +++ b/examples/advection/gas_dynamics_initials.py @@ -0,0 +1,172 @@ +# Copyright (C) 2008 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 . + + + + +from __future__ import division +from __future__ import absolute_import +import numpy +import numpy.linalg as la +from six.moves import range +from grudge import sym + + + + +class UniformMachFlow: + def __init__(self, mach=0.1, p=1, rho=1, reynolds=100, + gamma=1.4, prandtl=0.72, char_length=1, spec_gas_const=287.1, + angle_of_attack=None, direction=None, gaussian_pulse_at=None, + pulse_magnitude=0.1): + """ + :param direction: is a vector indicating the direction of the + flow. Only one of angle_of_attack and direction may be + specified. Only the direction, not the magnitude, of + direction is taken into account. + + :param angle_of_attack: if not None, specifies the angle of + the flow along the Y axis, where the flow is + directed along the X axis. + """ + if angle_of_attack is not None and direction is not None: + raise ValueError("Only one of angle_of_attack and " + "direction may be specified.") + + if angle_of_attack is None and direction is None: + angle_of_attack = 0 + + if direction is not None: + self.direction = direction/la.norm(direction) + else: + self.direction = None + + self.mach = mach + self.p = p + self.rho = rho + + self.gamma = gamma + self.prandtl = prandtl + self.reynolds = reynolds + self.length = char_length + self.spec_gas_const = spec_gas_const + + self.angle_of_attack = angle_of_attack + + self.gaussian_pulse_at = gaussian_pulse_at + self.pulse_magnitude = pulse_magnitude + + self.c = (self.gamma * p / rho)**0.5 + u = self.velocity = mach * self.c + self.e = p / (self.gamma - 1) + rho / 2 * u**2 + + if numpy.isinf(self.reynolds): + self.mu = 0 + else: + self.mu = u * self.length * rho / self.reynolds + + def direction_vector(self, dimensions): + # this must be done here because dimensions is not known above + if self.direction is None: + assert self.angle_of_attack is not None + direction = numpy.zeros(dimensions, dtype=numpy.float64) + direction[0] = numpy.cos( + self.angle_of_attack / 180. * numpy.pi) + direction[1] = numpy.sin( + self.angle_of_attack / 180. * numpy.pi) + return direction + else: + return self.direction + + def __call__(self, t, x_vec): + ones = numpy.ones_like(x_vec[0]) + rho_field = ones*self.rho + + if self.gaussian_pulse_at is not None: + rel_to_pulse = [x_vec[i] - self.gaussian_pulse_at[i] + for i in range(len(x_vec))] + rho_field += self.pulse_magnitude * self.rho * numpy.exp( + - sum(rtp_i**2 for rtp_i in rel_to_pulse)/2) + + direction = self.direction_vector(x_vec.shape[0]) + + from grudge.tools import make_obj_array + u_field = make_obj_array([ones*self.velocity*dir_i + for dir_i in direction]) + + from grudge.tools import join_fields + return join_fields(rho_field, self.e*ones, self.rho*u_field) + + def volume_interpolant(self, t, discr): + return discr.convert_volume( + self(t, discr.nodes.T), + kind=discr.compute_kind, + dtype=discr.default_scalar_type) + + def boundary_interpolant(self, t, discr, tag): + return discr.convert_boundary( + self(t, discr.get_boundary(tag).nodes.T), + tag=tag, kind=discr.compute_kind, + dtype=discr.default_scalar_type) + + + + + + + +class Vortex: + def __init__(self): + self.beta = 5 + self.gamma = 1.4 + self.center = numpy.array([5, 0]) + self.velocity = numpy.array([0.3, 0.3]) + self.final_time = 0.5 + + self.mu = 0 + self.prandtl = 0.72 + self.spec_gas_const = 287.1 + + def __call__(self): + x_vec = sym.nodes(2) + t = sym.var("t", sym.DD_SCALAR) + vortex_loc = self.center + t*self.velocity + + # coordinates relative to vortex center + x_rel = x_vec[0] - vortex_loc[0] + y_rel = x_vec[1] - vortex_loc[1] + + # Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159 + # also JSH/TW Nodal DG Methods, p. 209 + + from math import pi + r = sym.sqrt(x_rel**2+y_rel**2) + expterm = self.beta*sym.exp(1-r**2) + u = self.velocity[0] - expterm*y_rel/(2*pi) + v = self.velocity[1] + expterm*x_rel/(2*pi) + rho = (1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1)) + p = rho**self.gamma + + e = p/(self.gamma-1) + rho/2*(u**2+v**2) + + from pytools.obj_array import join_fields + return join_fields(rho, e, rho*u, rho*v) + + def volume_interpolant(self): + return self() + + + + diff --git a/examples/advection/new.py b/examples/advection/new.py new file mode 100644 index 00000000..6a8566e0 --- /dev/null +++ b/examples/advection/new.py @@ -0,0 +1,138 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2017 Bogdan Enache" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + +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 +from pytools.obj_array import join_fields + + +def main(write_output=True, order=4): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + dim = 2 + + resolution = 15 + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh(a=(0, -5), b=(10, 5), + n=(resolution, resolution), order=order) + + h = 1/resolution + + + source_center = np.array([3, -2]) + source_center2 = np.array([6, 2]) + source_width = 0.5 + + sym_x = sym.nodes(2) + sym_source_center_dist = sym_x - source_center + sym_source_center_dist2 = sym_x - source_center2 + + + def f_gaussian(x): + return join_fields(sym.exp( + -np.dot(sym_source_center_dist, sym_source_center_dist) + / source_width**2) ,sym.exp( + -np.dot(sym_source_center_dist2, sym_source_center_dist2) + / source_width**2) ) + + from gas_dynamics_initials import Vortex + flow = Vortex() + bund = sym.interp("vol", sym.BTAG_ALL)(flow.volume_interpolant()[0:2]) + + + + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) + + def sym_operator(): + u = sym.make_sym_array("u", 2) + def mflux(pair): + normal = sym.normal(pair.dd, dim) + return pair.avg*np.dot(normal, [1,1]) + + def flux(pair): + return sym.interp(pair.dd, "all_faces")( + mflux(pair)) + stiff_t = sym.stiffness_t(dim) + return sym.InverseMassOperator()( + np.dot(stiff_t, [u,u]) + - sym.FaceMassOperator()( + flux(sym.int_tpair(u)) + + flux(sym.bv_tpair(sym.BTAG_ALL, u, bund)) + + ) + ) + + + bound_op = bind(discr, sym_operator()) + + #u = bind(discr, f_gaussian(sym.nodes(dim)))(queue, t=0) + u = bind(discr, flow.volume_interpolant()[0:2])(queue, t=0) + + def rhs(t, u): + return bound_op(queue, t=t, u=u) + + final_time = 1 + dt = 0.001 + print(dt) + + from grudge.shortcuts import set_up_rk4 + dt_stepper = set_up_rk4("u", dt, u, rhs) + + from grudge.shortcuts import make_visualizer + vis = make_visualizer(discr, vis_order=2*order) + + step = 0 + vis.write_vtk_file("fld-%04d.vtu" % step, + #[("u", event.state_component[0])]) + [("u", u[0]), ("u2", u[1])] + ) + + for event in dt_stepper.run(t_end=final_time): + if isinstance(event, dt_stepper.StateComputed): + + step += 1 + print(step) + + vis.write_vtk_file("fld-%04d.vtu" % step, + #[("u", event.state_component[0])]) + [("u", event.state_component[0]), ("u2", event.state_component[1])] + ) + + +if __name__ == "__main__": + main() diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index 00e4f0c8..7cf72489 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -56,9 +56,9 @@ def main(write_output=True, order=4): sym_x = sym.nodes(2) - advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2 + advec_v = sym_x/sym_x - flux_type = "upwind" + flux_type = "central" source_center = np.array([0.1, 0.1]) source_width = 0.05 @@ -109,6 +109,8 @@ def main(write_output=True, order=4): dt = dt_factor * h/order**2 nsteps = (final_time // dt) + 1 dt = final_time/nsteps + 1e-15 + dt = 0.001 + print(dt) from grudge.shortcuts import set_up_rk4 dt_stepper = set_up_rk4("u", dt, u, rhs) @@ -122,7 +124,7 @@ def main(write_output=True, order=4): if isinstance(event, dt_stepper.StateComputed): step += 1 - if step % 30 == 0: + if step % 2 == 0: print(step) vis.write_vtk_file("fld-%04d.vtu" % step, diff --git a/examples/gas_dynamics/euler/vortex2.py b/examples/gas_dynamics/euler/vortex2.py index 4a4871cd..6a4d425f 100644 --- a/examples/gas_dynamics/euler/vortex2.py +++ b/examples/gas_dynamics/euler/vortex2.py @@ -38,8 +38,8 @@ def main(write_output=True, visualize=False): from meshmode.mesh.generation import generate_regular_rect_mesh - for order in [3, 4, 5]: - mesh = generate_regular_rect_mesh(a=(0,-5), b=(10,5), n=(6,6), order=order) + for order in [5]: + mesh = generate_regular_rect_mesh(a=(0,-5), b=(10,5), n=(5,5), order=order) from gas_dynamics_initials import Vortex flow = Vortex() @@ -67,7 +67,7 @@ def main(write_output=True, visualize=False): bound_flow_vol = bind(discr, flow.volume_interpolant()) - fields = bound_flow_vol(queue, t=1) + fields = bound_flow_vol(queue, t=0) print(fields) euler_ex = bind(discr, op.sym_operator()) @@ -75,7 +75,7 @@ def main(write_output=True, visualize=False): max_eigval = [0] def rhs(t, q): #ode_rhs, speed = euler_ex(queue, t=t, q=q) - ting = euler_ex(queue, t=t, q=q) + return euler_ex(queue, t=t, q=q) ode_rhs = ting[:-1] speed = ting[-1] #pu.db @@ -93,8 +93,8 @@ def main(write_output=True, visualize=False): #stepper = SSP3TimeStepper( #vector_primitive_factory=discr.get_vector_primitive_factory()) from grudge.shortcuts import set_up_rk4 - dt = 0.001 - final_t = 0.1 + dt = 0.01 + final_t = 2 dt_stepper = set_up_rk4("q", dt, fields, rhs) if visualize: diff --git a/examples/gas_dynamics/gas_dynamics_initials.py b/examples/gas_dynamics/gas_dynamics_initials.py index 9d24aa7f..588764f1 100644 --- a/examples/gas_dynamics/gas_dynamics_initials.py +++ b/examples/gas_dynamics/gas_dynamics_initials.py @@ -132,7 +132,7 @@ class Vortex: self.beta = 5 self.gamma = 1.4 self.center = numpy.array([5, 0]) - self.velocity = numpy.array([1, 0]) + self.velocity = numpy.array([0.3, 0.3]) self.final_time = 0.5 self.mu = 0 @@ -142,13 +142,13 @@ class Vortex: def __call__(self): x_vec = sym.nodes(2) t = sym.var("t", sym.DD_SCALAR) + from pytools.obj_array import join_fields vortex_loc = self.center + t*self.velocity # coordinates relative to vortex center x_rel = x_vec[0] - vortex_loc[0] y_rel = x_vec[1] - vortex_loc[1] - - # Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159 +# Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159 # also JSH/TW Nodal DG Methods, p. 209 from math import pi @@ -161,7 +161,6 @@ class Vortex: e = p/(self.gamma-1) + rho/2*(u**2+v**2) - from pytools.obj_array import join_fields return join_fields(rho, e, rho*u, rho*v) def volume_interpolant(self): diff --git a/grudge/execution.py b/grudge/execution.py index b7e24d48..5c7228a6 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -114,7 +114,7 @@ class ExecutionMapper(mappers.Evaluator, cached_name += "abs" else: cached_name += expr.function.name - print(cached_name + ":") + #print(cached_name + ":") @memoize_in(self.bound_op, cached_name) def knl(): @@ -129,7 +129,7 @@ class ExecutionMapper(mappers.Evaluator, assert len(args) == 1 evt, (out,) = knl()(self.queue, a=args[0]) - print(out) + #print(out) return out diff --git a/grudge/models/advection.py b/grudge/models/advection.py index 091b9870..d496760c 100644 --- a/grudge/models/advection.py +++ b/grudge/models/advection.py @@ -165,36 +165,34 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator): raise ValueError("invalid flux type") def sym_operator(self): + def mflux(pair): + normal = sym.normal(pair.dd, self.ambient_dim) + return pair.avg*np.dot(normal, [1,1]) + + def flux(pair): + return sym.interp(pair.dd, "all_faces")( + mflux(pair)) u = sym.var("u") + stiff_t = sym.stiffness_t(self.ambient_dim) + return sym.InverseMassOperator()( + np.dot(stiff_t, [u,u]) + - sym.FaceMassOperator()( + flux(sym.int_tpair(u)) + )) # boundary conditions ------------------------------------------------- bc_in = self.inflow_u - all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag) boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag) - def flux(pair): - return sym.interp(pair.dd, all_faces_dd)( - self.flux(pair)) quad_dd = sym.DOFDesc("vol", self.quad_tag) to_quad = sym.interp("vol", quad_dd) - stiff_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol") quad_v = to_quad(self.v) quad_u = to_quad(u) - return sym.InverseMassOperator()( - (stiff_t[0](quad_u * quad_v[0]) + stiff_t[1](quad_u * quad_v[1])) - - sym.FaceMassOperator(all_faces_dd, "vol")( - flux(sym.int_tpair(u, self.quad_tag)) - + flux(sym.bv_tpair(boundary_dd, u, bc_in)) - - # FIXME: Add back support for inflow/outflow tags - #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in)) - #+ flux(sym.bv_tpair(self.outflow_tag, u, bc_out)) - )) # }}} # vim: foldmethod=marker diff --git a/grudge/models/euler2.py b/grudge/models/euler2.py index c7fb2224..5fb2663d 100644 --- a/grudge/models/euler2.py +++ b/grudge/models/euler2.py @@ -87,6 +87,7 @@ class EulerOperator(HyperbolicOperator): ): self.dimensions = dimensions self.equation_of_state = equation_of_state + self.bc_inflow = sym.interp("vol",sym.BTAG_ALL)(bc_inflow) def rho(self, q): @@ -109,6 +110,18 @@ class EulerOperator(HyperbolicOperator): def flux(self, q): from pytools import delta + return [ # one entry for each flux direction + join_fields( + # flux rho + self.rho(q), + + # flux E + 0, + + # flux rho_u + 0, 0 + ) + for i in range(self.dimensions)] return [ # one entry for each flux direction join_fields( @@ -137,11 +150,16 @@ class EulerOperator(HyperbolicOperator): fluxes_ph = flux from pymbolic.primitives import Max - #penalty = Max((wave_speed_ph.int,wave_speed_ph.ext))*(state_ph.ext-state_ph.int) + #penalty = Max((wave_speed_ph.int,wave_speed_ph.ext))*(normal.ext*state_ph.ext+normal.int*state_ph.int) penalty = 0 - num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) - - penalty) + return state.avg*np.dot([1,1],normal) + #dotty = np.dot(normal, flux) + return (sum(n_i*(sym.int_tpair(f_i).avg) for n_i, f_i in zip(normal, fluxes_ph))) + + #num_flux = 0.5*(sum(n_i*(sym.int_tpair(f_i).int+sym.int_tpair(f_i).ext) for n_i, f_i in zip(normal, fluxes_ph)) + num_flux = (sum(n_i*(sym.int_tpair(f_i).avg) for n_i, f_i in zip(normal, fluxes_ph)) + + penalty) return num_flux def characteristic_velocity(self,state): @@ -158,6 +176,21 @@ class EulerOperator(HyperbolicOperator): def sym_operator(self): q = sym.make_sym_array("q", self.dimensions+2) + def mflux(pair): + normal = sym.normal(pair.dd, self.dimensions) + return pair.avg*np.dot(normal, [1,1]) + + def flux(pair): + return sym.interp(pair.dd, "all_faces")( + mflux(pair)) + stiff_t = sym.stiffness_t(self.dimensions) + return sym.InverseMassOperator()( + np.dot(stiff_t, [q,q]) + - sym.FaceMassOperator()( + flux(sym.int_tpair(q)) + + flux(sym.bv_tpair(sym.BTAG_ALL, q,self.bc_inflow )) + )) + def to_all_faces(pair, dd): @@ -167,9 +200,9 @@ class EulerOperator(HyperbolicOperator): volq_flux = self.flux(q) stiff_t = sym.stiffness_t(self.dimensions) ans = sym.InverseMassOperator()( - np.dot(stiff_t, volq_flux) + np.dot(stiff_t, [q,q]) - sym.FaceMassOperator()( - to_all_faces(self.lax_flux(sym.int_tpair(speed), sym.int_tpair(q), sym.int_tpair(self.flux(q)) ), sym.int_tpair(speed).dd) + to_all_faces(self.lax_flux(sym.int_tpair(speed), sym.int_tpair(q), self.flux(q) ), sym.int_tpair(q).dd) #+ flux(sym.bv_tpair(sym.BTAG_ALL, q, bc_in)) ) -- GitLab From 0d5cae97aa0a2e729aa92a8316d3eee6b4ffe6ff Mon Sep 17 00:00:00 2001 From: Bogdan Enache Date: Mon, 27 Aug 2018 19:51:00 -0500 Subject: [PATCH 15/15] Added some comments, made most recent commit be working --- examples/gas_dynamics/euler/vortex2.py | 23 +------------------ .../gas_dynamics/gas_dynamics_initials.py | 2 +- examples/maxwell/cavities.py | 4 ++-- grudge/models/gas_dynamics/__init__.py | 10 ++++---- 4 files changed, 8 insertions(+), 31 deletions(-) diff --git a/examples/gas_dynamics/euler/vortex2.py b/examples/gas_dynamics/euler/vortex2.py index 6a4d425f..81f19c29 100644 --- a/examples/gas_dynamics/euler/vortex2.py +++ b/examples/gas_dynamics/euler/vortex2.py @@ -15,7 +15,6 @@ - from __future__ import division from __future__ import absolute_import from __future__ import print_function @@ -41,7 +40,7 @@ def main(write_output=True, visualize=False): for order in [5]: mesh = generate_regular_rect_mesh(a=(0,-5), b=(10,5), n=(5,5), order=order) from gas_dynamics_initials import Vortex - flow = Vortex() + flow = Vortex() # NOT ACTUALLY A VORTEX, constant velocity from grudge.models.euler2 import EulerOperator from grudge.models.gas_dynamics import ( @@ -83,15 +82,6 @@ def main(write_output=True, visualize=False): return ode_rhs rhs(0, fields) - - # limiter ------------------------------------------------------------ - #from grudge.models.gas_dynamics import SlopeLimiter1NEuler - #limiter = SlopeLimiter1NEuler(discr, flow.gamma, 2, op) - - #from grudge.timestep.runge_kutta import SSP3TimeStepper - #stepper = SSP3TimeStepper(limiter=limiter) - #stepper = SSP3TimeStepper( - #vector_primitive_factory=discr.get_vector_primitive_factory()) from grudge.shortcuts import set_up_rk4 dt = 0.01 final_t = 2 @@ -160,16 +150,5 @@ def main(write_output=True, visualize=False): # after order loop assert eoc_rec.estimate_order_of_convergence()[0,1] > 6 - - - if __name__ == "__main__": main(visualize=True) - - - -# entry points for py.test ---------------------------------------------------- -from pytools.test import mark_test -@mark_test.long -def test_euler_vortex(): - main(write_output=False) diff --git a/examples/gas_dynamics/gas_dynamics_initials.py b/examples/gas_dynamics/gas_dynamics_initials.py index 588764f1..98cf671b 100644 --- a/examples/gas_dynamics/gas_dynamics_initials.py +++ b/examples/gas_dynamics/gas_dynamics_initials.py @@ -127,7 +127,7 @@ class UniformMachFlow: -class Vortex: +class Vortex: # NOT ACTUALLY A VORTEX, constant velocity def __init__(self): self.beta = 5 self.gamma = 1.4 diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py index 85a7c3de..a8763e84 100644 --- a/examples/maxwell/cavities.py +++ b/examples/maxwell/cavities.py @@ -28,7 +28,7 @@ THE SOFTWARE. import numpy as np import pyopencl as cl from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, Discretization +from grudge import sym, bind, DGDiscretizationWithBoundaries from grudge.models.em import get_rectangular_cavity_mode @@ -43,7 +43,7 @@ def main(dims, write_output=True, order=4): b=(1.0,)*dims, n=(5,)*dims) - discr = Discretization(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) if 0: epsilon0 = 8.8541878176e-12 # C**2 / (N m**2) diff --git a/grudge/models/gas_dynamics/__init__.py b/grudge/models/gas_dynamics/__init__.py index 34783833..4ff7e83a 100644 --- a/grudge/models/gas_dynamics/__init__.py +++ b/grudge/models/gas_dynamics/__init__.py @@ -557,13 +557,11 @@ class GasDynamicsOperator(TimeDependentOperator): # {{{ operator template --------------------------------------------------- def make_lax_friedrichs_flux(self, wave_speed, state, fluxes, bdry_tags_states_and_fluxes,strong): - #from hedge.flux import FluxVectorPlaceholder, flux_max def lff(fvph): n = len(state) d = len(fluxes) - normal = sym.normal(fvph.dd, d) #TODO: Also contradiction between this and some dd errors at the next todo - #fvph = FluxVectorPlaceholder(len(state)*(1+d)+1) + normal = sym.normal(fvph.dd, d) wave_speed_ph = fvph[0] state_ph = fvph[1:1+n] @@ -575,10 +573,10 @@ class GasDynamicsOperator(TimeDependentOperator): penalty = 0 if not strong: - num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) # TODO dd error between flux and norm here + num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) - penalty) else: - num_flux = 0.5*(sum(n_i*(f_i.int-f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) # TODO (or here) + num_flux = 0.5*(sum(n_i*(f_i.int-f_i.ext) for n_i, f_i in zip(normal, fluxes_ph)) + penalty) return num_flux @@ -589,7 +587,7 @@ class GasDynamicsOperator(TimeDependentOperator): int_operand = join_fields(wave_speed, state, *fluxes) - return fluxy(sym.int_tpair(int_operand, self.quad_tag)) + #return fluxy(sym.int_tpair(int_operand, self.quad_tag)) return (fluxy(sym.int_tpair(int_operand, self.quad_tag)) + sum( fluxy(sym.bv_tpair(sym.DOFDesc(sym.DTAG_BOUNDARY(tag), self.quad_tag),int_operand,join_fields(1, bdry_state, *bdry_fluxes) )) for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes)) -- GitLab