From 5e528b24628110e25125328d23dccfd5e97cfee7 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 11:43:50 -0600 Subject: [PATCH 01/25] Add projection/interpolation operators --- grudge/interpolation.py | 180 ++++++++++++++++++++++++++++++++++++++++ grudge/projection.py | 70 +++++++++++++++- 2 files changed, 249 insertions(+), 1 deletion(-) diff --git a/grudge/interpolation.py b/grudge/interpolation.py index 61bdf1a13..754f76abf 100644 --- a/grudge/interpolation.py +++ b/grudge/interpolation.py @@ -32,8 +32,22 @@ """ +import numpy as np + +from arraycontext import ( + ArrayContext, + map_array_container, + thaw +) +from functools import partial +from meshmode.transform_metadata import FirstAxisIsElementsTag + from grudge.discretization import DiscretizationCollection +from meshmode.dof_array import DOFArray + +from pytools import keyed_memoize_in + # FIXME: Should revamp interp and make clear distinctions # between projection and interpolations. @@ -46,3 +60,169 @@ def interp(dcoll: DiscretizationCollection, src, tgt, vec): from grudge.projection import project return project(dcoll, src, tgt, vec) + + +# {{{ Interpolation matrices + +def volume_quadrature_interpolation_matrix( + actx: ArrayContext, base_element_group, vol_quad_element_group): + """todo. + """ + @keyed_memoize_in( + actx, volume_quadrature_interpolation_matrix, + lambda base_grp, vol_quad_grp: (base_grp.discretization_key(), + vol_quad_grp.discretization_key())) + def get_volume_vand(base_grp, vol_quad_grp): + from modepy import vandermonde + + basis = base_grp.basis_obj() + vdm_inv = np.linalg.inv(vandermonde(basis.functions, + base_grp.unit_nodes)) + vdm_q = vandermonde(basis.functions, vol_quad_grp.unit_nodes) @ vdm_inv + return actx.freeze(actx.from_numpy(vdm_q)) + + return get_volume_vand(base_element_group, vol_quad_element_group) + + +def surface_quadrature_interpolation_matrix( + actx: ArrayContext, + base_element_group, face_quad_element_group, dtype): + """todo. + """ + @keyed_memoize_in( + actx, surface_quadrature_interpolation_matrix, + lambda base_grp, face_quad_grp: (base_grp.discretization_key(), + face_quad_grp.discretization_key())) + def get_surface_vand(base_grp, face_quad_grp): + nfaces = base_grp.mesh_el_group.nfaces + assert face_quad_grp.nelements == nfaces * base_grp.nelements + + from modepy import vandermonde, faces_for_shape + + basis = base_grp.basis_obj() + vdm_inv = np.linalg.inv(vandermonde(basis.functions, + base_grp.unit_nodes)) + faces = faces_for_shape(base_grp.shape) + # NOTE: Assumes same quadrature rule on each face + face_quadrature = face_quad_grp.quadrature_rule() + + surface_nodes = faces[0].map_to_volume(face_quadrature.nodes) + for fidx in range(1, nfaces): + surface_nodes = np.append( + surface_nodes, + faces[fidx].map_to_volume(face_quadrature.nodes), + axis=1 + ) + vdm_f = vandermonde(basis.functions, surface_nodes) @ vdm_inv + return actx.freeze(actx.from_numpy(vdm_f)) + + return get_surface_vand(base_element_group, face_quad_element_group) + + +def volume_and_surface_interpolation_matrix( + actx: ArrayContext, + base_element_group, + vol_quad_element_group, + face_quad_element_group, dtype): + """todo. + """ + @keyed_memoize_in( + actx, volume_and_surface_interpolation_matrix, + lambda base_grp, vol_quad_grp, face_quad_grp: ( + base_grp.discretization_key(), + vol_quad_grp.discretization_key(), + face_quad_grp.discretization_key())) + def get_vol_surf_interpolation_matrix(base_grp, vol_quad_grp, face_quad_grp): + vq_mat = actx.to_numpy( + thaw( + volume_quadrature_interpolation_matrix( + actx, base_grp, vol_quad_grp + ), + actx + ) + ) + vf_mat = actx.to_numpy( + thaw( + surface_quadrature_interpolation_matrix( + actx, base_grp, face_quad_grp, dtype + ), + actx + ) + ) + return actx.freeze(actx.from_numpy(np.block([[vq_mat], [vf_mat]]))) + + return get_vol_surf_interpolation_matrix( + base_element_group, vol_quad_element_group, face_quad_element_group + ) + +# }}} + + +def volume_quadrature_interpolation(dcoll: DiscretizationCollection, dq, vec): + if not isinstance(vec, DOFArray): + return map_array_container( + partial(volume_quadrature_interpolation, dcoll, dq), vec + ) + + actx = vec.array_context + discr = dcoll.discr_from_dd("vol") + quad_discr = dcoll.discr_from_dd(dq) + + return DOFArray( + actx, + data=tuple( + actx.einsum("ij,ej->ei", + volume_quadrature_interpolation_matrix( + actx, + base_element_group=bgrp, + vol_quad_element_group=qgrp + ), + vec_i, + arg_names=("Vq_mat", "vec"), + tagged=(FirstAxisIsElementsTag(),)) + + for bgrp, qgrp, vec_i in zip(discr.groups, quad_discr.groups, vec) + ) + ) + + +def volume_and_surface_quadrature_interpolation( + dcoll: DiscretizationCollection, dq, df, vec): + """todo. + """ + if not isinstance(vec, DOFArray): + return map_array_container( + partial(volume_and_surface_quadrature_interpolation, + dcoll, dq, df), vec + ) + + actx = vec.array_context + dtype = vec.entry_dtype + discr = dcoll.discr_from_dd("vol") + quad_volm_discr = dcoll.discr_from_dd(dq) + quad_face_discr = dcoll.discr_from_dd(df) + + return DOFArray( + actx, + data=tuple( + actx.einsum("ij,ej->ei", + volume_and_surface_interpolation_matrix( + actx, + base_element_group=bgrp, + vol_quad_element_group=qvgrp, + face_quad_element_group=qfgrp, + dtype=dtype + ), + vec_i, + arg_names=("Vh_mat", "vec"), + tagged=(FirstAxisIsElementsTag(),)) + + for bgrp, qvgrp, qfgrp, vec_i in zip( + discr.groups, + quad_volm_discr.groups, + quad_face_discr.groups, vec) + ) + ) + + +# vim: foldmethod=marker diff --git a/grudge/projection.py b/grudge/projection.py index fdf951555..0380f245f 100644 --- a/grudge/projection.py +++ b/grudge/projection.py @@ -32,15 +32,20 @@ """ +import numpy as np + from functools import partial -from arraycontext import map_array_container +from arraycontext import ArrayContext, map_array_container from arraycontext.container import ArrayOrContainerT from grudge.discretization import DiscretizationCollection from grudge.dof_desc import as_dofdesc from meshmode.dof_array import DOFArray +from meshmode.transform_metadata import FirstAxisIsElementsTag + +from pytools import keyed_memoize_in from numbers import Number @@ -70,3 +75,66 @@ def project( ) return dcoll.connection_from_dds(src, tgt)(vec) + + +# {{{ Projection matrices + +def volume_quadrature_l2_projection_matrix( + actx: ArrayContext, base_element_group, vol_quad_element_group): + """todo. + """ + @keyed_memoize_in( + actx, volume_quadrature_l2_projection_matrix, + lambda base_grp, vol_quad_grp: (base_grp.discretization_key(), + vol_quad_grp.discretization_key())) + def get_ref_l2_proj_mat(base_grp, vol_quad_grp): + from grudge.interpolation import volume_quadrature_interpolation_matrix + + vdm_q = actx.to_numpy( + volume_quadrature_interpolation_matrix( + actx, base_grp, vol_quad_grp + ) + ) + weights = vol_quad_grp.quadrature_rule().weights + inv_mass_mat = np.linalg.inv(weights * vdm_q.T @ vdm_q) + return actx.freeze(actx.from_numpy(inv_mass_mat @ (vdm_q.T * weights))) + + return get_ref_l2_proj_mat(base_element_group, vol_quad_element_group) + +# }}} + + +def volume_quadrature_project(dcoll: DiscretizationCollection, dd_q, vec): + """todo. + """ + if not isinstance(vec, DOFArray): + return map_array_container( + partial(volume_quadrature_project, dcoll, dd_q), vec + ) + + actx = vec.array_context + discr = dcoll.discr_from_dd("vol") + quad_discr = dcoll.discr_from_dd(dd_q) + + return DOFArray( + actx, + data=tuple( + actx.einsum("ij,ej->ei", + volume_quadrature_l2_projection_matrix( + actx, + base_element_group=bgrp, + vol_quad_element_group=qgrp + ), + vec_i, + arg_names=("Pq_mat", "vec"), + tagged=(FirstAxisIsElementsTag(),)) + + for bgrp, qgrp, vec_i in zip( + discr.groups, + quad_discr.groups, + vec + ) + ) + ) + +# vim: foldmethod=marker From 8cf491273916c772763e7d778b540c5061ccdf66 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 11:44:03 -0600 Subject: [PATCH 02/25] Add unit test for SBP properties --- test/test_sbp_ops.py | 175 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 175 insertions(+) create mode 100644 test/test_sbp_ops.py diff --git a/test/test_sbp_ops.py b/test/test_sbp_ops.py new file mode 100644 index 000000000..05e8b8ee6 --- /dev/null +++ b/test/test_sbp_ops.py @@ -0,0 +1,175 @@ +__copyright__ = "Copyright (C) 2021 University of Illinois Board of Trustees" + +__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 + +from grudge import DiscretizationCollection +from grudge.dof_desc import DOFDesc, DISCR_TAG_BASE, DISCR_TAG_QUAD + +import pytest + +from grudge.array_context import PytestPyOpenCLArrayContextFactory +from arraycontext import pytest_generate_tests_for_array_contexts +pytest_generate_tests = pytest_generate_tests_for_array_contexts( + [PytestPyOpenCLArrayContextFactory]) + +import logging + +logger = logging.getLogger(__name__) + + +@pytest.mark.parametrize("dim", [1, 2, 3]) +@pytest.mark.parametrize("order", [1, 2, 3, 4]) +def test_reference_element_sbp_operators(actx_factory, dim, order): + actx = actx_factory() + + from meshmode.mesh.generation import generate_regular_rect_mesh + + nel_1d = 5 + box_ll = -5.0 + box_ur = 5.0 + mesh = generate_regular_rect_mesh( + a=(box_ll,)*dim, + b=(box_ur,)*dim, + nelements_per_axis=(nel_1d,)*dim) + + from meshmode.discretization.poly_element import \ + default_simplex_group_factory, QuadratureSimplexGroupFactory + + dcoll = DiscretizationCollection( + actx, mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: default_simplex_group_factory(dim, order), + DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order) + } + ) + + dd_q = DOFDesc("vol", DISCR_TAG_QUAD) + dd_f = DOFDesc("all_faces", DISCR_TAG_QUAD) + + volm_discr = dcoll.discr_from_dd("vol") + quad_discr = dcoll.discr_from_dd(dd_q) + quad_face_discr = dcoll.discr_from_dd(dd_f) + dtype = volm_discr.zeros(actx).entry_dtype + + from meshmode.discretization.poly_element import diff_matrices + from modepy import faces_for_shape, face_normal + from grudge.projection import volume_quadrature_l2_projection_matrix + from grudge.interpolation import ( + volume_quadrature_interpolation_matrix, + surface_quadrature_interpolation_matrix + ) + + for vgrp, qgrp, qfgrp in zip(volm_discr.groups, + quad_discr.groups, + quad_face_discr.groups): + nq_vol = qgrp.nunit_dofs + nq_per_face = qfgrp.nunit_dofs + nfaces = vgrp.shape.nfaces + nq_faces = nfaces * nq_per_face + nq_total = nq_vol + nq_faces + + # {{{ Volume operators + + weights = qgrp.quadrature_rule().weights + vdm_q = actx.to_numpy( + volume_quadrature_interpolation_matrix(actx, vgrp, qgrp) + ) + mass_mat = weights * vdm_q.T @ vdm_q + p_mat = actx.to_numpy( + volume_quadrature_l2_projection_matrix(actx, vgrp, qgrp) + ) + + # }}} + + # Checks Pq @ Vq = Minv @ Vq.T @ W @ Vq = I + assert np.allclose(p_mat @ vdm_q, + np.identity(len(mass_mat)), rtol=1e-15) + + # {{{ Surface operators + + faces = faces_for_shape(vgrp.shape) + # NOTE: assumes same quadrature rule on all faces + face_weights = np.tile(qfgrp.quadrature_rule().weights, nfaces) + face_normals = [face_normal(face) for face in faces] + e = np.ones(shape=(nq_per_face,)) + nrstj = [np.concatenate([np.sign(nhat[idx])*e + for nhat in face_normals]) + for idx in range(vgrp.dim)] + b_mats = [np.diag(face_weights*nrstj[d]) for d in range(vgrp.dim)] + vf_mat = actx.to_numpy( + surface_quadrature_interpolation_matrix( + actx, + base_element_group=vgrp, + face_quad_element_group=qfgrp, + dtype=dtype + ) + ) + + # }}} + + # {{{ Hybridized (volume + surface) operators + + q_mats = [p_mat.T @ mass_mat @ diff_mat @ p_mat + for diff_mat in diff_matrices(vgrp)] + e_mat = vf_mat @ p_mat + qtilde_mats = 0.5 * np.asarray( + [np.block([[q_mats[d] - q_mats[d].T, e_mat.T @ b_mats[d]], + [-b_mats[d] @ e_mat, b_mats[d]]]) + for d in range(dcoll.dim)] + ) + + # }}} + + ones = np.ones(shape=(nq_total,)) + zeros = np.zeros(shape=(nq_total,)) + for idx in range(dim): + # Checks the generalized SBP property: + # Qi + Qi.T = E.T @ Bi @ E + # c.f. Lemma 1. https://arxiv.org/pdf/1708.01243.pdf + assert np.allclose(q_mats[idx] + q_mats[idx].T, + e_mat.T @ b_mats[idx] @ e_mat, rtol=1e-15) + + # Checks the SBP-like property for the skew hybridized operator + # Qitilde + Qitilde.T = [0 0; 0 Bi] + # c.f. Theorem 1 and Lemma 1. https://arxiv.org/pdf/1902.01828.pdf + residual = qtilde_mats[idx] + qtilde_mats[idx].T + residual[nq_vol:nq_vol+nq_faces, nq_vol:nq_vol+nq_faces] -= b_mats[idx] + assert np.allclose(residual, np.zeros(residual.shape), rtol=1e-15) + + # Checks quadrature condition for: Qiskew @ ones = zeros + # Qiskew + Qiskew.T = [0 0; 0 Bi] + # c.f. Lemma 2. https://arxiv.org/pdf/1902.01828.pdf + assert np.allclose(np.dot(qtilde_mats[idx], ones), + zeros, rtol=1e-15) + + +# You can test individual routines by typing +# $ python test_grudge.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) From fca7f4678914996e9fe89fb1e3a3485495c058d0 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 14:05:34 -0600 Subject: [PATCH 03/25] Use SingleGridWorkBalancingPytatoArrayContext --- grudge/array_context.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/grudge/array_context.py b/grudge/array_context.py index 0f85a62bb..2a624f863 100644 --- a/grudge/array_context.py +++ b/grudge/array_context.py @@ -28,7 +28,10 @@ from meshmode.array_context import ( PyOpenCLArrayContext as _PyOpenCLArrayContextBase, - PytatoPyOpenCLArrayContext as _PytatoPyOpenCLArrayContextBase, + # FIXME: PytatoPyOpenCLArrayContext is literally unusable when + # using flux-differencing + # TODO: Get SingleGridWorkBalancingPytatoArrayContext merged into main + SingleGridWorkBalancingPytatoArrayContext as _PytatoPyOpenCLArrayContextBase, ) from arraycontext.pytest import ( _PytestPyOpenCLArrayContextFactoryWithClass, From 7499e9f82bc06d791382a3d1ee2a31c0b1363a2a Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 14:06:14 -0600 Subject: [PATCH 04/25] Add entropy stable Euler model and flux-differencing --- grudge/flux_differencing.py | 272 ++++++++++++++++++++++++++++++++++++ grudge/models/euler.py | 252 ++++++++++++++++++++++++++++++++- 2 files changed, 523 insertions(+), 1 deletion(-) create mode 100644 grudge/flux_differencing.py diff --git a/grudge/flux_differencing.py b/grudge/flux_differencing.py new file mode 100644 index 000000000..01036c42c --- /dev/null +++ b/grudge/flux_differencing.py @@ -0,0 +1,272 @@ +__copyright__ = """ +Copyright (C) 2021 University of Illinois Board of Trustees +""" + +__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. +""" + + +from arraycontext import ( + ArrayContext, + map_array_container, + freeze +) +from arraycontext.container import ArrayOrContainerT +from typing import Callable +from functools import partial +from meshmode.transform_metadata import FirstAxisIsElementsTag + +from grudge.discretization import DiscretizationCollection + +from meshmode.dof_array import DOFArray + +from pytools import memoize_in, keyed_memoize_in + +import numpy as np + + +def skew_symmetric_hybridized_sbp_operators( + actx: ArrayContext, + base_element_group, + vol_quad_element_group, + face_quad_element_group, dtype): + """todo. + """ + @keyed_memoize_in( + actx, skew_symmetric_hybridized_sbp_operators, + lambda base_grp, quad_vol_grp, face_quad_grp: ( + base_grp.discretization_key(), + quad_vol_grp.discretization_key(), + face_quad_grp.discretization_key() + ) + ) + def get_skew_symetric_hybridized_diff_mats( + base_grp, quad_vol_grp, face_quad_grp): + from meshmode.discretization.poly_element import diff_matrices + from modepy import faces_for_shape, face_normal + from grudge.projection import volume_quadrature_l2_projection_matrix + from grudge.interpolation import ( + volume_quadrature_interpolation_matrix, + surface_quadrature_interpolation_matrix + ) + + # {{{ Volume operators + + weights = quad_vol_grp.quadrature_rule().weights + vdm_q = actx.to_numpy( + volume_quadrature_interpolation_matrix(actx, base_grp, quad_vol_grp) + ) + mass_mat = weights * vdm_q.T @ vdm_q + p_mat = actx.to_numpy( + volume_quadrature_l2_projection_matrix(actx, base_grp, quad_vol_grp) + ) + + # }}} + + # {{{ Surface operators + + faces = faces_for_shape(base_grp.shape) + nfaces = len(faces) + # NOTE: assumes same quadrature rule on all faces + face_weights = np.tile(face_quad_grp.quadrature_rule().weights, nfaces) + face_normals = [face_normal(face) for face in faces] + nnods_per_face = face_quad_grp.nunit_dofs + e = np.ones(shape=(nnods_per_face,)) + nrstj = [ + # nsrtJ = nhat * Jhatf, where nhat is the reference normal + # and Jhatf is the Jacobian det. of the transformation from + # the face of the reference element to the reference face. + np.concatenate([np.sign(nhat[idx])*e for nhat in face_normals]) + for idx in range(base_grp.dim) + ] + b_mats = [np.diag(face_weights*nrstj[d]) for d in range(base_grp.dim)] + vf_mat = actx.to_numpy( + surface_quadrature_interpolation_matrix( + actx, + base_element_group=base_grp, + face_quad_element_group=face_quad_grp, + dtype=dtype + ) + ) + zero_mat = np.zeros((nfaces*nnods_per_face, nfaces*nnods_per_face), + dtype=dtype) + + # }}} + + # {{{ Hybridized (volume + surface) operators + + q_mats = [p_mat.T @ mass_mat @ diff_mat @ p_mat + for diff_mat in diff_matrices(base_grp)] + e_mat = vf_mat @ p_mat + q_skew_hybridized = np.asarray( + [ + np.block( + [[q_mats[d] - q_mats[d].T, e_mat.T @ b_mats[d]], + [-b_mats[d] @ e_mat, zero_mat]] + ) for d in range(base_grp.dim) + ], + order="C" + ) + + # }}} + + return actx.freeze(actx.from_numpy(q_skew_hybridized)) + + return get_skew_symetric_hybridized_diff_mats( + base_element_group, + vol_quad_element_group, + face_quad_element_group + ) + + +def _single_axis_hybridized_sbp_derivative_kernel( + dcoll, dq, df, xyz_axis, flux_matrix): + if not isinstance(flux_matrix, DOFArray): + return map_array_container( + partial(_single_axis_hybridized_sbp_derivative_kernel, + dcoll, dq, df, xyz_axis), + flux_matrix + ) + + from grudge.geometry import \ + area_element, inverse_surface_metric_derivative + from grudge.interpolation import ( + volume_and_surface_interpolation_matrix, + volume_and_surface_quadrature_interpolation + ) + + # FIXME: This is kinda meh + def inverse_jac_matrix(): + @memoize_in(dcoll, (_single_axis_hybridized_sbp_derivative_kernel, dq, df)) + def _inv_surf_metric_deriv(): + mat = actx.np.stack( + [ + actx.np.stack( + [ + volume_and_surface_quadrature_interpolation( + dcoll, dq, df, + area_element(actx, dcoll) + * inverse_surface_metric_derivative( + actx, dcoll, + rst_ax, xyz_ax + ) + ) for rst_ax in range(dcoll.dim) + ] + ) for xyz_ax in range(dcoll.ambient_dim) + ] + ) + return freeze(mat, actx) + return _inv_surf_metric_deriv() + + actx = flux_matrix.array_context + discr = dcoll.discr_from_dd("vol") + quad_volm_discr = dcoll.discr_from_dd(dq) + quad_face_discr = dcoll.discr_from_dd(df) + + # FIXME: This assumes affine-geometry + return DOFArray( + actx, + data=tuple( + # r for rst axis + actx.einsum("ik,rej,rij,eij->ek", + volume_and_surface_interpolation_matrix( + actx, + base_element_group=bgrp, + vol_quad_element_group=qvgrp, + face_quad_element_group=qafgrp, + dtype=fmat_i.dtype + ), + ijm_i[xyz_axis], + skew_symmetric_hybridized_sbp_operators( + actx, + bgrp, + qvgrp, + qafgrp, + fmat_i.dtype + ), + fmat_i, + arg_names=("Vh_mat_t", "inv_jac_t", "ref_Q_mat", "F_mat"), + tagged=(FirstAxisIsElementsTag(),)) + + for bgrp, qvgrp, qafgrp, fmat_i, ijm_i in zip( + discr.groups, + quad_volm_discr.groups, + quad_face_discr.groups, + flux_matrix, + inverse_jac_matrix() + ) + ) + ) + + +def _apply_skew_symmetric_hybrid_diff_operator( + dcoll: DiscretizationCollection, dq, df, flux_matrices): + # flux matrices for each component are a vector of matrices (as DOFArrays) + # for each spatial dimension + if not (isinstance(flux_matrices, np.ndarray) and flux_matrices.dtype == "O"): + return map_array_container( + partial(_apply_skew_symmetric_hybrid_diff_operator, dcoll, dq, df), + flux_matrices + ) + + def _sbp_hybrid_diff_helper(diff_func, mats): + if not isinstance(mats, np.ndarray): + raise TypeError("argument must be an object array") + assert mats.dtype == object + return sum(diff_func(i, mat_i) for i, mat_i in enumerate(mats)) + + return _sbp_hybrid_diff_helper( + lambda i, flux_mat_i: _single_axis_hybridized_sbp_derivative_kernel( + dcoll, dq, df, i, flux_mat_i), + flux_matrices + ) + + +def volume_flux_differencing( + dcoll: DiscretizationCollection, + flux: Callable[[ArrayOrContainerT, ArrayOrContainerT], ArrayOrContainerT], + dq, df, vec: ArrayOrContainerT) -> ArrayOrContainerT: + """todo. + """ + def _reshape_to_numpy(shape, ary): + if not isinstance(ary, DOFArray): + return map_array_container(partial(_reshape_to_numpy, shape), ary) + + actx = ary.array_context + return DOFArray( + actx, + data=tuple( + subary.reshape(grp.nelements, *shape) + for grp, subary in zip( + dcoll.discr_from_dd("vol").groups, + ary + ) + ) + ) + + return _apply_skew_symmetric_hybrid_diff_operator( + dcoll, + dq, df, + flux(_reshape_to_numpy((1, -1), vec), + _reshape_to_numpy((-1, 1), vec)) + ) + + +# vim: foldmethod=marker diff --git a/grudge/models/euler.py b/grudge/models/euler.py index e3cc0b6b5..db8a83dc2 100644 --- a/grudge/models/euler.py +++ b/grudge/models/euler.py @@ -142,7 +142,65 @@ def conservative_to_primitive_vars(cv_state: ConservedEulerField, gamma=1.4): u = rho_u / rho p = (gamma - 1) * (rho_e - 0.5 * sum(rho_u * u)) - return rho, u, p + return (rho, u, p) + + +def conservative_to_entropy_vars(cv_state, gamma=1.4): + """Converts from conserved variables (density, momentum, total energy) + into entropy variables. + + :arg cv_state: A :class:`EulerContainer` containing the conserved + variables. + :arg gamma: The isentropic expansion factor for a single-species gas + (default set to 1.4). + :returns: A :class:`EulerContainer` containing the entropy variables. + """ + actx = cv_state.array_context + rho, u, p = conservative_to_primitive_vars(cv_state, gamma=gamma) + + u_square = sum(v ** 2 for v in u) + s = actx.np.log(p) - gamma*actx.np.log(rho) + rho_p = rho / p + + return ConservedEulerField( + mass=((gamma - s)/(gamma - 1)) - 0.5 * rho_p * u_square, + energy=-rho_p, + momentum=rho_p * u + ) + + +def entropy_to_conservative_vars(ev_state, gamma=1.4): + """Converts from entropy variables into conserved variables + (density, momentum, total energy). + + :arg ev_state: A :class:`EulerContainer` containing the entropy + variables. + :arg gamma: The isentropic expansion factor for a single-species gas + (default set to 1.4). + :returns: A :class:`EulerContainer` containing the conserved variables. + """ + actx = ev_state.array_context + # See Hughes, Franca, Mallet (1986) A new finite element + # formulation for CFD: (DOI: 10.1016/0045-7825(86)90127-1) + inv_gamma_minus_one = 1/(gamma - 1) + + # Convert to entropy `-rho * s` used by Hughes, France, Mallet (1986) + ev_state = ev_state * (gamma - 1) + v1 = ev_state.mass + v2t4 = ev_state.momentum + v5 = ev_state.energy + + v_square = sum(v**2 for v in v2t4) + s = gamma - v1 + v_square/(2*v5) + rho_iota = ( + ((gamma - 1) / (-v5)**gamma)**(inv_gamma_minus_one) + ) * actx.np.exp(-s * inv_gamma_minus_one) + + return ConservedEulerField( + mass=-rho_iota * v5, + energy=rho_iota * (1 - v_square/(2*v5)), + momentum=rho_iota * v2t4 + ) def compute_wavespeed(cv_state: ConservedEulerField, gamma=1.4): @@ -373,4 +431,196 @@ def interp_to_quad_surf(u): # }}} +# {{{ Entropy stable Euler operator + +def flux_chandrashekar(dcoll: DiscretizationCollection, q_ll, q_rr, gamma=1.4): + """Two-point volume flux based on the entropy conserving + and kinetic energy preserving two-point flux in: + + - Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite + Volume Schemes for Compressible Euler and Navier-Stokes Equations + [DOI](https://doi.org/10.4208/cicp.170712.010313a) + + :args q_ll: an array container for the "left" state. + :args q_rr: an array container for the "right" state. + :arg gamma: The isentropic expansion factor for a single-species gas + (default set to 1.4). + """ + dim = dcoll.dim + actx = q_ll.array_context + + def ln_mean(x: DOFArray, y: DOFArray, epsilon=1e-4): + f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) + return actx.np.where( + actx.np.less(f2, epsilon), + (x + y) / (2 + f2*2/3 + f2*f2*2/5 + f2*f2*f2*2/7), + (y - x) / actx.np.log(y / x) + ) + + rho_ll, u_ll, p_ll = conservative_to_primitive_vars(q_ll, gamma=gamma) + rho_rr, u_rr, p_rr = conservative_to_primitive_vars(q_rr, gamma=gamma) + + beta_ll = 0.5 * rho_ll / p_ll + beta_rr = 0.5 * rho_rr / p_rr + specific_kin_ll = 0.5 * sum(v**2 for v in u_ll) + specific_kin_rr = 0.5 * sum(v**2 for v in u_rr) + + rho_avg = 0.5 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) + beta_mean = ln_mean(beta_ll, beta_rr) + beta_avg = 0.5 * (beta_ll + beta_rr) + u_avg = 0.5 * (u_ll + u_rr) + p_mean = 0.5 * rho_avg / beta_avg + + velocity_square_avg = specific_kin_ll + specific_kin_rr + + mass_flux = rho_mean * u_avg + momentum_flux = np.outer(mass_flux, u_avg) + np.eye(dim) * p_mean + energy_flux = ( + mass_flux * 0.5 * (1/(gamma - 1)/beta_mean - velocity_square_avg) + + np.dot(momentum_flux, u_avg) + ) + + return EulerContainer(mass=mass_flux, + energy=energy_flux, + momentum=momentum_flux) + + +def entropy_stable_numerical_flux_chandrashekar( + dcoll: DiscretizationCollection, tpair: TracePair, + gamma=1.4, lf_stabilization=False): + """Entropy stable numerical flux based on the entropy conserving + and kinetic energy preserving two-point flux in: + + - Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite + Volume Schemes for Compressible Euler and Navier-Stokes Equations + [DOI](https://doi.org/10.4208/cicp.170712.010313a) + """ + dd_intfaces = tpair.dd + dd_allfaces = dd_intfaces.with_dtag("all_faces") + q_ll = tpair.int + q_rr = tpair.ext + actx = q_ll.array_context + + num_flux = flux_chandrashekar(dcoll, q_ll, q_rr, gamma=gamma) + normal = thaw(dcoll.normal(dd_intfaces), actx) + + if lf_stabilization: + from arraycontext import outer + + rho_ll, u_ll, p_ll = conservative_to_primitive_vars(q_ll, gamma=gamma) + rho_rr, u_rr, p_rr = conservative_to_primitive_vars(q_rr, gamma=gamma) + + def compute_wavespeed(rho, u, p): + return ( + actx.np.sqrt(np.dot(u, u)) + actx.np.sqrt(gamma * (p / rho)) + ) + + # Compute jump penalization parameter + lam = actx.np.maximum(compute_wavespeed(rho_ll, u_ll, p_ll), + compute_wavespeed(rho_rr, u_rr, p_rr)) + num_flux -= lam*outer(tpair.diff, normal)/2 + + return op.project(dcoll, dd_intfaces, dd_allfaces, num_flux @ normal) + + +class EntropyStableEulerOperator(EulerOperator): + + def operator(self, t, q): + from grudge.projection import volume_quadrature_project + from grudge.interpolation import \ + volume_and_surface_quadrature_interpolation + + gamma = self.gamma + qtag = self.qtag + dq = DOFDesc("vol", qtag) + df = DOFDesc("all_faces", qtag) + + dcoll = self.dcoll + + # Interpolate cv_state to vol quad grid: u_q = V_q u + q_quad = op.project(dcoll, "vol", dq, q) + + # Convert to projected entropy variables: v_q = V_h P_q v(u_q) + entropy_vars = conservative_to_entropy_vars(q_quad, gamma=gamma) + proj_entropy_vars = volume_quadrature_project(dcoll, dq, entropy_vars) + + # Compute conserved state in terms of the (interpolated) + # projected entropy variables on the quad grid + qtilde_allquad = \ + entropy_to_conservative_vars( + # Interpolate projected entropy variables to + # volume + surface quadrature grids + volume_and_surface_quadrature_interpolation( + dcoll, dq, df, proj_entropy_vars + ), + gamma=gamma) + + # Compute volume derivatives using flux differencing + from functools import partial + from grudge.flux_differencing import volume_flux_differencing + + flux_diff = volume_flux_differencing( + dcoll, + partial(flux_chandrashekar, dcoll, gamma=gamma), + dq, df, + qtilde_allquad + ) + + def entropy_tpair(tpair): + dd_intfaces = tpair.dd + dd_intfaces_quad = dd_intfaces.with_discr_tag(qtag) + # Interpolate entropy variables to the surface quadrature grid + vtilde_tpair = \ + op.project(dcoll, dd_intfaces, dd_intfaces_quad, tpair) + return TracePair( + dd_intfaces_quad, + # Convert interior and exterior states to conserved variables + interior=entropy_to_conservative_vars( + vtilde_tpair.int, gamma=gamma + ), + exterior=entropy_to_conservative_vars( + vtilde_tpair.ext, gamma=gamma + ) + ) + + # Computing interface numerical fluxes + interface_fluxes = ( + sum( + entropy_stable_numerical_flux_chandrashekar( + dcoll, + entropy_tpair(tpair), + gamma=gamma, + lf_stabilization=self.lf_stabilization + # NOTE: Trace pairs consist of the projected entropy variables + # which will be converted to conserved variables in + # *entropy_tpair* + ) for tpair in op.interior_trace_pairs(dcoll, proj_entropy_vars) + ) + ) + + # Compute boundary fluxes + if self.bdry_conditions is not None: + bc_fluxes = sum( + entropy_stable_numerical_flux_chandrashekar( + dcoll, + self.bdry_conditions[btag].boundary_tpair( + dcoll, + as_dofdesc(btag).with_discr_tag(qtag), + q, + t=t + ), + gamma=gamma, + lf_stabilization=self.lf_stabilization + ) for btag in self.bdry_conditions + ) + interface_fluxes = interface_fluxes + bc_fluxes + + return op.inverse_mass( + dcoll, + -flux_diff - op.face_mass(dcoll, df, interface_fluxes) + ) + +# }}} + # vim: foldmethod=marker From cb555a5422c0470fb5c034f3386a527e0ac473d8 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 14:06:35 -0600 Subject: [PATCH 05/25] Test variable transformations in Euler module --- test/test_euler_module.py | 134 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 134 insertions(+) create mode 100644 test/test_euler_module.py diff --git a/test/test_euler_module.py b/test/test_euler_module.py new file mode 100644 index 000000000..79059ae9f --- /dev/null +++ b/test/test_euler_module.py @@ -0,0 +1,134 @@ +__copyright__ = "Copyright (C) 2021 University of Illinois Board of Trustees" + +__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 + +from arraycontext import thaw + +from grudge import DiscretizationCollection +from grudge.dof_desc import DISCR_TAG_BASE + +import grudge.op as op + +from pytools.obj_array import make_obj_array + +import pytest + +from grudge.array_context import PytestPyOpenCLArrayContextFactory +from arraycontext import pytest_generate_tests_for_array_contexts +pytest_generate_tests = pytest_generate_tests_for_array_contexts( + [PytestPyOpenCLArrayContextFactory]) + +import logging + +logger = logging.getLogger(__name__) + + +def test_variable_transformations(actx_factory): + from grudge.models.euler import ( + primitive_to_conservative_vars, + conservative_to_primitive_vars, + entropy_to_conservative_vars, + conservative_to_entropy_vars + ) + + actx = actx_factory() + gamma = 1.4 # Adiabatic expansion factor for single-gas Euler model + + def vortex(x_vec, t=0): + mach = 0.5 # Mach number + _x0 = 5 + epsilon = 1 # vortex strength + x, y = x_vec + actx = x.array_context + + fxyt = 1 - (((x - _x0) - t)**2 + y**2) + expterm = actx.np.exp(fxyt/2) + + u = 1 - (epsilon*y/(2*np.pi))*expterm + v = ((epsilon*(x - _x0) - t)/(2*np.pi))*expterm + + velocity = make_obj_array([u, v]) + mass = ( + 1 - ((epsilon**2*(gamma - 1)*mach**2)/(8*np.pi**2))*actx.np.exp(fxyt) + ) ** (1 / (gamma - 1)) + p = (mass**gamma)/(gamma*mach**2) + + return primitive_to_conservative_vars((mass, velocity, p), gamma=gamma) + + from meshmode.mesh.generation import generate_regular_rect_mesh + + dim = 2 + res = 5 + mesh = generate_regular_rect_mesh( + a=(0, -5), + b=(20, 5), + nelements_per_axis=(2*res, res), + periodic=(True, True)) + + from meshmode.discretization.poly_element import \ + default_simplex_group_factory + + order = 3 + dcoll = DiscretizationCollection( + actx, mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: default_simplex_group_factory(dim, order) + } + ) + + # Fields in conserved variables + fields = vortex(thaw(dcoll.nodes(), actx)) + + # Map back and forth between primitive and conserved vars + fields_prim = conservative_to_primitive_vars(fields, gamma=gamma) + prim_fields_to_cons = primitive_to_conservative_vars(fields_prim, gamma=gamma) + + assert op.norm( + dcoll, prim_fields_to_cons.mass - fields.mass, np.inf) < 1e-13 + assert op.norm( + dcoll, prim_fields_to_cons.energy - fields.energy, np.inf) < 1e-13 + assert op.norm( + dcoll, prim_fields_to_cons.momentum - fields.momentum, np.inf) < 1e-13 + + # Map back and forth between entropy and conserved vars + fields_ev = conservative_to_entropy_vars(fields, gamma=gamma) + ev_fields_to_cons = entropy_to_conservative_vars(fields_ev, gamma=gamma) + + assert op.norm( + dcoll, ev_fields_to_cons.mass - fields.mass, np.inf) < 1e-13 + assert op.norm( + dcoll, ev_fields_to_cons.energy - fields.energy, np.inf) < 1e-13 + assert op.norm( + dcoll, ev_fields_to_cons.momentum - fields.momentum, np.inf) < 1e-13 + + +# You can test individual routines by typing +# $ python test_grudge.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) From b0eb50a77005c28663b99e0a5187fb9caf325bdc Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 14:06:51 -0600 Subject: [PATCH 06/25] Set up ESDG in vortex example --- examples/euler/vortex.py | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/examples/euler/vortex.py b/examples/euler/vortex.py index 0e4cf7d76..f712d0082 100644 --- a/examples/euler/vortex.py +++ b/examples/euler/vortex.py @@ -31,7 +31,8 @@ from grudge.array_context import PytatoPyOpenCLArrayContext, PyOpenCLArrayContext from grudge.models.euler import ( vortex_initial_condition, - EulerOperator + EulerOperator, + EntropyStableEulerOperator ) from grudge.shortcuts import rk4_step @@ -42,6 +43,7 @@ def run_vortex(actx, order=3, resolution=8, final_time=5, + esdg=False, overintegration=False, flux_type="central", visualize=False): @@ -52,11 +54,12 @@ def run_vortex(actx, order=3, resolution=8, final_time=5, order: %s\n final_time: %s\n resolution: %s\n + entropy stable: %s\n overintegration: %s\n flux_type: %s\n visualize: %s """, - order, final_time, resolution, + order, final_time, resolution, esdg, overintegration, flux_type, visualize ) @@ -78,7 +81,14 @@ def run_vortex(actx, order=3, resolution=8, final_time=5, from meshmode.discretization.poly_element import \ default_simplex_group_factory, QuadratureSimplexGroupFactory - exp_name = f"fld-vortex-N{order}-K{resolution}-{flux_type}" + if esdg: + case = "esdg-vortex" + operator_cls = EntropyStableEulerOperator + else: + case = "vortex" + operator_cls = EulerOperator + + exp_name = f"fld-{case}-N{order}-K{resolution}-{flux_type}" if overintegration: exp_name += "-overintegrated" @@ -99,7 +109,7 @@ def run_vortex(actx, order=3, resolution=8, final_time=5, # {{{ Euler operator - euler_operator = EulerOperator( + euler_operator = operator_cls( dcoll, flux_type=flux_type, gamma=gamma, @@ -156,6 +166,7 @@ def rhs(t, q): def main(ctx_factory, order=3, final_time=5, resolution=8, + esdg=False, overintegration=False, lf_stabilization=False, visualize=False, @@ -175,6 +186,12 @@ def main(ctx_factory, order=3, final_time=5, resolution=8, force_device_scalars=True, ) + if not actx.supports_nonscalar_broadcasting and esdg is True: + raise RuntimeError( + "Cannot use ESDG with an array context that cannot perform " + "nonscalar broadcasting. Run with --lazy instead." + ) + if lf_stabilization: flux_type = "lf" else: @@ -185,6 +202,7 @@ def main(ctx_factory, order=3, final_time=5, resolution=8, order=order, resolution=resolution, final_time=final_time, + esdg=esdg, overintegration=overintegration, flux_type=flux_type, visualize=visualize @@ -198,6 +216,8 @@ def main(ctx_factory, order=3, final_time=5, resolution=8, parser.add_argument("--order", default=3, type=int) parser.add_argument("--tfinal", default=0.015, type=float) parser.add_argument("--resolution", default=8, type=int) + parser.add_argument("--esdg", action="store_true", + help="use entropy stable dg") parser.add_argument("--oi", action="store_true", help="use overintegration") parser.add_argument("--lf", action="store_true", @@ -213,6 +233,7 @@ def main(ctx_factory, order=3, final_time=5, resolution=8, order=args.order, final_time=args.tfinal, resolution=args.resolution, + esdg=args.esdg, overintegration=args.oi, lf_stabilization=args.lf, visualize=args.visualize, From 956d3058ca67031d4788c7e8eb7f1baa1173433f Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 14:22:36 -0600 Subject: [PATCH 07/25] Update requirements.txt to point to the development branches --- requirements.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index abf6d39fa..ccd240b99 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,13 +4,13 @@ git+https://github.com/inducer/pytools.git#egg=pytools git+https://github.com/inducer/pymbolic.git#egg=pymbolic git+https://github.com/inducer/islpy.git#egg=islpy git+https://github.com/inducer/pyopencl.git#egg=pyopencl -git+https://github.com/inducer/loopy.git#egg=loopy +git+https://github.com/kaushikcfd/loopy.git@pytato-array-context-transforms#egg=loopy git+https://github.com/inducer/dagrt.git#egg=dagrt git+https://github.com/inducer/leap.git#egg=leap git+https://github.com/inducer/meshpy.git#egg=meshpy git+https://github.com/inducer/modepy.git#egg=modepy git+https://github.com/inducer/arraycontext.git#egg=arraycontext -git+https://github.com/inducer/meshmode.git#egg=meshmode +git+https://github.com/kaushikcfd/meshmode.git@pytato-array-context-transforms#egg=meshmode git+https://github.com/inducer/pyvisfile.git#egg=pyvisfile git+https://github.com/inducer/pymetis.git#egg=pymetis git+https://github.com/illinois-ceesd/logpyle.git#egg=logpyle From b510e27c111c02abc5098ce10a16c521230c2c78 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 14:39:50 -0600 Subject: [PATCH 08/25] Add esdg option for pulse problem --- examples/euler/acoustic_pulse.py | 40 +++++++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/examples/euler/acoustic_pulse.py b/examples/euler/acoustic_pulse.py index 2672dda12..2af79d16c 100644 --- a/examples/euler/acoustic_pulse.py +++ b/examples/euler/acoustic_pulse.py @@ -36,6 +36,7 @@ from grudge.models.euler import ( ConservedEulerField, EulerOperator, + EntropyStableEulerOperator, InviscidWallBC ) from grudge.shortcuts import rk4_step @@ -112,9 +113,24 @@ def run_acoustic_pulse(actx, order=3, final_time=1, resolution=16, + esdg=False, overintegration=False, visualize=False): + logger.info( + """ + Acoustic pulse parameters:\n + order: %s\n + final_time: %s\n + resolution: %s\n + entropy stable: %s\n + overintegration: %s\n + visualize: %s + """, + order, final_time, resolution, esdg, + overintegration, visualize + ) + # eos-related parameters gamma = 1.4 @@ -136,7 +152,15 @@ def run_acoustic_pulse(actx, (default_simplex_group_factory, QuadratureSimplexGroupFactory) - exp_name = f"fld-acoustic-pulse-N{order}-K{resolution}" + if esdg: + case = "esdg-pulse" + operator_cls = EntropyStableEulerOperator + else: + case = "pulse" + operator_cls = EulerOperator + + exp_name = f"fld-{case}-N{order}-K{resolution}" + if overintegration: exp_name += "-overintegrated" quad_tag = DISCR_TAG_QUAD @@ -156,7 +180,7 @@ def run_acoustic_pulse(actx, # {{{ Euler operator - euler_operator = EulerOperator( + euler_operator = operator_cls( dcoll, bdry_conditions={BTAG_ALL: InviscidWallBC()}, flux_type="lf", @@ -213,7 +237,7 @@ def rhs(t, q): def main(ctx_factory, order=3, final_time=1, resolution=16, - overintegration=False, visualize=False, lazy=False): + esdg=False, overintegration=False, visualize=False, lazy=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) @@ -229,10 +253,17 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, force_device_scalars=True, ) + if not actx.supports_nonscalar_broadcasting and esdg is True: + raise RuntimeError( + "Cannot use ESDG with an array context that cannot perform " + "nonscalar broadcasting. Run with --lazy instead." + ) + run_acoustic_pulse( actx, order=order, resolution=resolution, + esdg=esdg, overintegration=overintegration, final_time=final_time, visualize=visualize @@ -246,6 +277,8 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, parser.add_argument("--order", default=3, type=int) parser.add_argument("--tfinal", default=0.1, type=float) parser.add_argument("--resolution", default=16, type=int) + parser.add_argument("--esdg", action="store_true", + help="use entropy stable dg") parser.add_argument("--oi", action="store_true", help="use overintegration") parser.add_argument("--visualize", action="store_true", @@ -259,6 +292,7 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, order=args.order, final_time=args.tfinal, resolution=args.resolution, + esdg=args.esdg, overintegration=args.oi, visualize=args.visualize, lazy=args.lazy) From cf7eeadd209b7919aaf147c4d56aae225d2de0cc Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 15:57:04 -0600 Subject: [PATCH 09/25] Add lsrk methods to shortcuts --- grudge/shortcuts.py | 118 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 118 insertions(+) diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py index 52c77097d..73f743a4e 100644 --- a/grudge/shortcuts.py +++ b/grudge/shortcuts.py @@ -22,6 +22,124 @@ THE SOFTWARE. """ +import numpy as np + +from dataclasses import dataclass + + +@dataclass(frozen=True) +class LSRKCoefficients: + """Dataclass which defines a given low-storage Runge-Kutta (LSRK) scheme. + The methods are determined by the provided `A`, `B` and `C` coefficient arrays. + """ + + A: np.ndarray + B: np.ndarray + C: np.ndarray + + +LSRK144NiegemannDiehlBuschCoefs = LSRKCoefficients( + A=np.array([ + 0., + -0.7188012108672410, + -0.7785331173421570, + -0.0053282796654044, + -0.8552979934029281, + -3.9564138245774565, + -1.5780575380587385, + -2.0837094552574054, + -0.7483334182761610, + -0.7032861106563359, + 0.0013917096117681, + -0.0932075369637460, + -0.9514200470875948, + -7.1151571693922548]), + B=np.array([ + 0.0367762454319673, + 0.3136296607553959, + 0.1531848691869027, + 0.0030097086818182, + 0.3326293790646110, + 0.2440251405350864, + 0.3718879239592277, + 0.6204126221582444, + 0.1524043173028741, + 0.0760894927419266, + 0.0077604214040978, + 0.0024647284755382, + 0.0780348340049386, + 5.5059777270269628]), + C=np.array([ + 0., + 0.0367762454319673, + 0.1249685262725025, + 0.2446177702277698, + 0.2476149531070420, + 0.2969311120382472, + 0.3978149645802642, + 0.5270854589440328, + 0.6981269994175695, + 0.8190890835352128, + 0.8527059887098624, + 0.8604711817462826, + 0.8627060376969976, + 0.8734213127600976])) + + +LSRK54CarpenterKennedyCoefs = LSRKCoefficients( + A=np.array([ + 0., + -567301805773/1357537059087, + -2404267990393/2016746695238, + -3550918686646/2091501179385, + -1275806237668/842570457699]), + B=np.array([ + 1432997174477/9575080441755, + 5161836677717/13612068292357, + 1720146321549/2090206949498, + 3134564353537/4481467310338, + 2277821191437/14882151754819]), + C=np.array([ + 0., + 1432997174477/9575080441755, + 2526269341429/6820363962896, + 2006345519317/3224310063776, + 2802321613138/2924317926251])) + + +EulerCoefs = LSRKCoefficients( + A=np.array([0.]), + B=np.array([1.]), + C=np.array([0.])) + + +def lsrk_step(coefs, state, t, dt, rhs): + """Take one step using a low-storage Runge-Kutta method.""" + k = 0.0 * state + for i in range(len(coefs.A)): + k = coefs.A[i]*k + dt*rhs(t + coefs.C[i]*dt, state) + state += coefs.B[i]*k + return state + + +def lsrk144_step(state, t, dt, rhs): + """Take one step using an explicit 14-stage, 4th-order, LSRK method. + + This method is derived by Niegemann, Diehl, and Busch (2012), with + an optimal stability region for advection-dominated flows. + """ + return lsrk_step(LSRK144NiegemannDiehlBuschCoefs, state, t, dt, rhs) + + +def lsrk54_step(state, t, dt, rhs): + """Take one step using an explicit 5-stage, 4th-order, LSRK method.""" + return lsrk_step(LSRK54CarpenterKennedyCoefs, state, t, dt, rhs) + + +def euler_step(state, t, dt, rhs): + """Take one step using the explicit, 1st-order accurate, Euler method.""" + return lsrk_step(EulerCoefs, state, t, dt, rhs) + def rk4_step(y, t, h, f): k1 = f(t, y) From 63b1ef0222f63fb4e3f7b7a4cf9f56de9d453e4b Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Dec 2021 15:57:35 -0600 Subject: [PATCH 10/25] Readd 1-D Sod shock problem --- examples/euler/sod.py | 271 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 271 insertions(+) create mode 100644 examples/euler/sod.py diff --git a/examples/euler/sod.py b/examples/euler/sod.py new file mode 100644 index 000000000..dda36cf71 --- /dev/null +++ b/examples/euler/sod.py @@ -0,0 +1,271 @@ +"""Minimal example of a grudge driver.""" + +__copyright__ = """ +Copyright (C) 2021 University of Illinois Board of Trustees +""" + +__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 pyopencl as cl +import pyopencl.tools as cl_tools + +from arraycontext import thaw, freeze +from meshmode.array_context import ( + SingleGridWorkBalancingPytatoArrayContext as PytatoPyOpenCLArrayContext +) +from grudge.models.euler import ( + EntropyStableEulerOperator, + EulerContainer, + PrescribedBC, + conservative_to_primitive_vars, + primitive_to_conservative_vars +) +from grudge.shortcuts import lsrk54_step + +from pytools.obj_array import make_obj_array + +import grudge.op as op + +import logging +logger = logging.getLogger(__name__) + + +# def sod_shock_initial_condition(nodes, t=0): +# gamma = 1.4 +# dim = len(nodes) +# x = nodes[0] +# actx = x.array_context +# zeros = 0*x + +# _x0 = 0.5 +# _rhoin = 1.0 +# _rhoout = 0.125 +# _pin = 1.0 +# _pout = 0.1 +# rhoin = zeros + _rhoin +# rhoout = zeros + _rhoout + +# x0 = zeros + _x0 +# sigma = 1e-13 +# weight = 0.5 * (1.0 - actx.np.tanh(1.0/sigma * (x - x0))) + +# rho = rhoout + (rhoin - rhoout)*weight +# p = _pout + (_pin - _pout)*weight +# u = make_obj_array([zeros for _ in range(dim)]) + +# return primitive_to_conservative_vars((rho, u, p), gamma=gamma) + +def sod_shock_initial_condition(nodes, t=0): + gamma = 1.4 + dim = len(nodes) + gmn1 = 1.0 / (gamma - 1.0) + x = nodes[0] + actx = x.array_context + zeros = 0*x + + _x0 = 0.5 + _rhoin = 1.0 + _rhoout = 0.125 + _pin = 1.0 + _pout = 0.1 + rhoin = zeros + _rhoin + rhoout = zeros + _rhoout + energyin = zeros + gmn1 * _pin + energyout = zeros + gmn1 * _pout + + x0 = zeros + _x0 + sigma = 1e-13 + weight = 0.5 * (1.0 - actx.np.tanh(1.0/sigma * (x - x0))) + + mass = rhoout + (rhoin - rhoout)*weight + energy = energyout + (energyin - energyout)*weight + momentum = make_obj_array([zeros for _ in range(dim)]) + + return EulerContainer(mass=mass, energy=energy, momentum=momentum) + + +def run_sod_shock_tube(actx, + order=4, + resolution=32, + final_time=0.2, + overintegration=False, + visualize=False): + + logger.info( + """ + Sod 1-D parameters:\n + order: %s\n + final_time: %s\n + resolution: %s\n + overintegration: %s\n + visualize: %s + """, + order, final_time, resolution, + overintegration, visualize + ) + + # eos-related parameters + gamma = 1.4 + + # {{{ discretization + + from meshmode.mesh.generation import generate_regular_rect_mesh + + dim = 1 + box_ll = 0.0 + box_ur = 1.0 + mesh = generate_regular_rect_mesh( + a=(box_ll,)*dim, + b=(box_ur,)*dim, + nelements_per_axis=(resolution,)*dim, + boundary_tag_to_face={ + "prescribed": ["+x", "-x"], + } + ) + + from grudge import DiscretizationCollection + from grudge.dof_desc import \ + DISCR_TAG_BASE, DISCR_TAG_QUAD, DTAG_BOUNDARY + from meshmode.discretization.poly_element import \ + (default_simplex_group_factory, + QuadratureSimplexGroupFactory) + + exp_name = f"fld-sod-1d-N{order}-K{resolution}" + + if overintegration: + exp_name += "-overintegrated" + quad_tag = DISCR_TAG_QUAD + else: + quad_tag = None + + dcoll = DiscretizationCollection( + actx, mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: default_simplex_group_factory(dim, order), + DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(order + 2) + } + ) + + # }}} + + # {{{ Euler operator + + dd_prescribe = DTAG_BOUNDARY("prescribed") + bcs = { + dd_prescribe: PrescribedBC(prescribed_state=sod_shock_initial_condition) + } + + euler_operator = EntropyStableEulerOperator( + dcoll, + bdry_conditions=bcs, + flux_type="lf", + gamma=gamma, + quadrature_tag=quad_tag + ) + + def rhs(t, q): + return euler_operator.operator(t, q) + + compiled_rhs = actx.compile(rhs) + + fields = sod_shock_initial_condition(thaw(dcoll.nodes(), actx)) + + from grudge.dt_utils import h_min_from_volume + + cfl = 0.01 + cn = 0.5*(order + 1)**2 + dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn + + logger.info("Timestep size: %g", dt) + + # }}} + + from grudge.shortcuts import make_visualizer + + vis = make_visualizer(dcoll) + + # {{{ time stepping + + step = 0 + t = 0.0 + while t < final_time: + if step % 10 == 0: + norm_q = actx.to_numpy(op.norm(dcoll, fields, 2)) + logger.info("[%04d] t = %.5f |q| = %.5e", step, t, norm_q) + if visualize: + rho, velocity, pressure = \ + conservative_to_primitive_vars(fields, gamma=gamma) + vis.write_vtk_file( + f"{exp_name}-{step:04d}.vtu", + [ + ("rho", rho), + ("energy", fields.energy), + ("momentum", fields.momentum), + ("velocity", velocity), + ("pressure", pressure) + ] + ) + assert norm_q < 10000 + + fields = thaw(freeze(fields, actx), actx) + fields = lsrk54_step(fields, t, dt, compiled_rhs) + t += dt + step += 1 + + # }}} + + +def main(ctx_factory, order=4, final_time=0.2, resolution=32, + overintegration=False, visualize=False): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PytatoPyOpenCLArrayContext( + queue, + allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)), + ) + + run_sod_shock_tube( + actx, order=order, + resolution=resolution, + final_time=final_time, + overintegration=overintegration, + visualize=visualize) + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser() + parser.add_argument("--order", default=4, type=int) + parser.add_argument("--tfinal", default=0.2, type=float) + parser.add_argument("--resolution", default=32, type=int) + parser.add_argument("--oi", action="store_true") + parser.add_argument("--visualize", action="store_true") + args = parser.parse_args() + + logging.basicConfig(level=logging.INFO) + main(cl.create_some_context, + order=args.order, + final_time=args.tfinal, + resolution=args.resolution, + overintegration=args.oi, + visualize=args.visualize) From eeadc0317ebe3f23fa0e8c8071955018402aaffb Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 9 Dec 2021 10:25:19 -0600 Subject: [PATCH 11/25] Taylor Green vortex experiment --- examples/euler/taylor_green.py | 240 +++++++++++++++++++++++++++++++++ 1 file changed, 240 insertions(+) create mode 100644 examples/euler/taylor_green.py diff --git a/examples/euler/taylor_green.py b/examples/euler/taylor_green.py new file mode 100644 index 000000000..24bf07044 --- /dev/null +++ b/examples/euler/taylor_green.py @@ -0,0 +1,240 @@ +"""Minimal example of a grudge driver.""" + +__copyright__ = """ +Copyright (C) 2021 University of Illinois Board of Trustees +""" + +__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 +import pyopencl.tools as cl_tools + +from arraycontext import thaw, freeze +from grudge.array_context import PytatoPyOpenCLArrayContext +from grudge.models.euler import ( + primitive_to_conservative_vars, + conservative_to_primitive_vars, + EntropyStableEulerOperator +) +from grudge.shortcuts import lsrk54_step, lsrk144_step + +from pytools.obj_array import make_obj_array + +import grudge.op as op + +import logging +logger = logging.getLogger(__name__) + + +def tg_vortex_initial_condition(x_vec, t=0): + # Parameters chosen from https://arxiv.org/pdf/1909.12546.pdf + M = 0.05 + L = 1. + gamma = 1.4 + v0 = 1. + p0 = 1. + rho0 = gamma * M ** 2 + + x, y, z = x_vec + actx = x.array_context + ones = 0*x + 1 + + p = p0 + rho0 * (v0 ** 2) / 16 * ( + actx.np.cos(2*x / L + actx.np.cos(2*y / L))) * actx.np.cos(2*z / L + 2) + u = v0 * actx.np.sin(x / L) * actx.np.cos(y / L) * actx.np.cos(z / L) + v = -v0 * actx.np.cos(x / L) * actx.np.sin(y / L) * actx.np.cos(z / L) + w = 0 * z + velocity = make_obj_array([u, v, w]) + rho = rho0 * ones + + return primitive_to_conservative_vars((rho, velocity, p), gamma=gamma) + + +def run_tg_vortex(actx, order=3, resolution=16, dtscaling=1, final_time=20, + dumpfreq=50, + overintegration=False, + timestepper="lsrk54", visualize=False): + logger.info( + """ + Taylor-Green vortex parameters:\n + order: %s, resolution: %s, dt scaling factor: %s, \n + final time: %s, \n + dumpfreq: %s, + timestepper: %s, visualize: %s + """, + order, resolution, dtscaling, + final_time, dumpfreq, timestepper, visualize + ) + + # eos-related parameters + gamma = 1.4 + + # {{{ discretization + + from meshmode.mesh.generation import generate_regular_rect_mesh + + dim = 3 + box_ll = -np.pi + box_ur = np.pi + mesh = generate_regular_rect_mesh( + a=(box_ll,)*dim, + b=(box_ur,)*dim, + nelements_per_axis=(resolution,)*dim, + periodic=(True,)*dim) + + from grudge import DiscretizationCollection + from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD + from meshmode.discretization.poly_element import \ + (default_simplex_group_factory, + QuadratureSimplexGroupFactory) + + exp_name = f"fld-esdg-taylorgreen-N{order}-K{resolution}" + + if overintegration: + exp_name += "-overintegrated" + quad_tag = DISCR_TAG_QUAD + else: + quad_tag = None + + dcoll = DiscretizationCollection( + actx, mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: default_simplex_group_factory(dim, order), + DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order) + } + ) + + # }}} + + # {{{ Euler operator + + euler_operator = EntropyStableEulerOperator( + dcoll, + flux_type="lf", + gamma=gamma, + quadrature_tag=quad_tag + ) + + def rhs(t, q): + return euler_operator.operator(t, q) + + compiled_rhs = actx.compile(rhs) + + fields = tg_vortex_initial_condition(thaw(dcoll.nodes(), actx)) + + if timestepper == "lsrk54": + stepper = lsrk54_step + else: + assert timestepper == "lsrk144" + stepper = lsrk144_step + + dt = dtscaling * actx.to_numpy( + euler_operator.estimate_rk4_timestep(actx, dcoll, state=fields)) + + logger.info("Timestep size: %g", dt) + + # }}} + + from grudge.shortcuts import make_visualizer + + vis = make_visualizer(dcoll) + + # {{{ time stepping + + step = 0 + t = 0.0 + while t < final_time: + if step % dumpfreq == 0: + norm_q = actx.to_numpy(op.norm(dcoll, fields, 2)) + logger.info("[%04d] t = %.5f |q| = %.5e", step, t, norm_q) + if visualize: + rho, velocity, pressure = \ + conservative_to_primitive_vars(fields, gamma=gamma) + vis.write_vtk_file( + f"fld-taylor-green-vortex-{step:04d}.vtu", + [ + ("rho", rho), + ("energy", fields.energy), + ("momentum", fields.momentum), + ("velocity", velocity), + ("pressure", pressure) + ] + ) + + fields = thaw(freeze(fields, actx), actx) + fields = stepper(fields, t, dt, compiled_rhs) + t += dt + step += 1 + + # }}} + + +def main(ctx_factory, order=3, final_time=20, + resolution=16, dtscaling=1, + dumpfreq=50, overintegration=False, + timestepper="lsrk54", visualize=False): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PytatoPyOpenCLArrayContext( + queue, + allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)) + ) + run_tg_vortex( + actx, + order=order, + resolution=resolution, + dtscaling=dtscaling, + final_time=final_time, + dumpfreq=dumpfreq, + overintegration=overintegration, + timestepper=timestepper, + visualize=visualize + ) + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser() + parser.add_argument("--order", default=3, type=int) + parser.add_argument("--tfinal", default=20., type=float) + parser.add_argument("--resolution", default=8, type=int) + parser.add_argument("--dumpfreq", default=50, type=int) + parser.add_argument("--dtscaling", default=1., type=float) + parser.add_argument("--oi", action="store_true") + parser.add_argument("--visualize", action="store_true") + parser.add_argument("--timestepper", default="lsrk54", + choices=['lsrk54', 'lsrk144']) + args = parser.parse_args() + + logging.basicConfig(level=logging.INFO) + main(cl.create_some_context, + order=args.order, + final_time=args.tfinal, + resolution=args.resolution, + dtscaling=args.dtscaling, + dumpfreq=args.dumpfreq, + overintegration=args.oi, + timestepper=args.timestepper, + visualize=args.visualize) From 8fe4f645f89cf5c4b1e1cf45584a116e5d75bc7d Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 9 Dec 2021 11:03:16 -0600 Subject: [PATCH 12/25] Clean up and document taylor green experiment --- examples/euler/taylor_green.py | 64 ++++++++++++++++++++++------------ 1 file changed, 42 insertions(+), 22 deletions(-) diff --git a/examples/euler/taylor_green.py b/examples/euler/taylor_green.py index 24bf07044..df80f644d 100644 --- a/examples/euler/taylor_green.py +++ b/examples/euler/taylor_green.py @@ -54,7 +54,7 @@ def tg_vortex_initial_condition(x_vec, t=0): gamma = 1.4 v0 = 1. p0 = 1. - rho0 = gamma * M ** 2 + rho0 = gamma * (M ** 2) x, y, z = x_vec actx = x.array_context @@ -71,20 +71,30 @@ def tg_vortex_initial_condition(x_vec, t=0): return primitive_to_conservative_vars((rho, velocity, p), gamma=gamma) -def run_tg_vortex(actx, order=3, resolution=16, dtscaling=1, final_time=20, +def run_tg_vortex(actx, + order=3, + resolution=16, + cfl=1.0, + final_time=20, dumpfreq=50, overintegration=False, - timestepper="lsrk54", visualize=False): + timestepper="lsrk54", + visualize=False): + logger.info( """ Taylor-Green vortex parameters:\n - order: %s, resolution: %s, dt scaling factor: %s, \n - final time: %s, \n - dumpfreq: %s, - timestepper: %s, visualize: %s + order: %s\n + final_time: %s\n + resolution: %s\n + cfl: %s\n + dumpfreq: %s\n + overintegration: %s\n + timestepper: %s\n + visualize: %s """, - order, resolution, dtscaling, - final_time, dumpfreq, timestepper, visualize + order, final_time, resolution, cfl, + dumpfreq, overintegration, timestepper, visualize ) # eos-related parameters @@ -109,7 +119,7 @@ def run_tg_vortex(actx, order=3, resolution=16, dtscaling=1, final_time=20, (default_simplex_group_factory, QuadratureSimplexGroupFactory) - exp_name = f"fld-esdg-taylorgreen-N{order}-K{resolution}" + exp_name = f"fld-esdg-taylorgreen-N{order}-K{resolution}-cfl{cfl}" if overintegration: exp_name += "-overintegrated" @@ -149,8 +159,10 @@ def rhs(t, q): assert timestepper == "lsrk144" stepper = lsrk144_step - dt = dtscaling * actx.to_numpy( - euler_operator.estimate_rk4_timestep(actx, dcoll, state=fields)) + from grudge.dt_utils import h_min_from_volume + + cn = 0.5*(order + 1)**2 + dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn logger.info("Timestep size: %g", dt) @@ -181,6 +193,7 @@ def rhs(t, q): ("pressure", pressure) ] ) + assert norm_q < 10000 fields = thaw(freeze(fields, actx), actx) fields = stepper(fields, t, dt, compiled_rhs) @@ -191,7 +204,7 @@ def rhs(t, q): def main(ctx_factory, order=3, final_time=20, - resolution=16, dtscaling=1, + resolution=16, cfl=1, dumpfreq=50, overintegration=False, timestepper="lsrk54", visualize=False): cl_ctx = ctx_factory() @@ -204,7 +217,7 @@ def main(ctx_factory, order=3, final_time=20, actx, order=order, resolution=resolution, - dtscaling=dtscaling, + cfl=cfl, final_time=final_time, dumpfreq=dumpfreq, overintegration=overintegration, @@ -218,14 +231,21 @@ def main(ctx_factory, order=3, final_time=20, parser = argparse.ArgumentParser() parser.add_argument("--order", default=3, type=int) - parser.add_argument("--tfinal", default=20., type=float) - parser.add_argument("--resolution", default=8, type=int) - parser.add_argument("--dumpfreq", default=50, type=int) - parser.add_argument("--dtscaling", default=1., type=float) - parser.add_argument("--oi", action="store_true") - parser.add_argument("--visualize", action="store_true") + parser.add_argument("--tfinal", default=20., type=float, + help="specify final time for the simulation") + parser.add_argument("--resolution", default=8, type=int, + help="resolution in each spatial direction") + parser.add_argument("--dumpfreq", default=50, type=int, + help="step frequency for writing output or logging info") + parser.add_argument("--cfl", default=1., type=float, + help="specify cfl value to use in dt computation") + parser.add_argument("--oi", action="store_true", + help="use overintegration") + parser.add_argument("--visualize", action="store_true", + help="write out vtk output") parser.add_argument("--timestepper", default="lsrk54", - choices=['lsrk54', 'lsrk144']) + choices=['lsrk54', 'lsrk144'], + help="specify timestepper method") args = parser.parse_args() logging.basicConfig(level=logging.INFO) @@ -233,7 +253,7 @@ def main(ctx_factory, order=3, final_time=20, order=args.order, final_time=args.tfinal, resolution=args.resolution, - dtscaling=args.dtscaling, + cfl=args.cfl, dumpfreq=args.dumpfreq, overintegration=args.oi, timestepper=args.timestepper, From 37acb3895e1fe651a08bbccca75f01e782f2ba7a Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Tue, 14 Dec 2021 10:46:53 -0600 Subject: [PATCH 13/25] Revamp helper function for flux-differencing --- grudge/flux_differencing.py | 61 ++++++++++++++++++++++--------------- 1 file changed, 37 insertions(+), 24 deletions(-) diff --git a/grudge/flux_differencing.py b/grudge/flux_differencing.py index 01036c42c..26a2abfdf 100644 --- a/grudge/flux_differencing.py +++ b/grudge/flux_differencing.py @@ -35,6 +35,8 @@ from grudge.discretization import DiscretizationCollection +import grudge.dof_desc as dof_desc + from meshmode.dof_array import DOFArray from pytools import memoize_in, keyed_memoize_in @@ -216,27 +218,38 @@ def _inv_surf_metric_deriv(): ) -def _apply_skew_symmetric_hybrid_diff_operator( - dcoll: DiscretizationCollection, dq, df, flux_matrices): - # flux matrices for each component are a vector of matrices (as DOFArrays) - # for each spatial dimension - if not (isinstance(flux_matrices, np.ndarray) and flux_matrices.dtype == "O"): +def _flux_differencing_helper(dcoll, diff_func, mats): + if not isinstance(mats, np.ndarray): + # mats is not an object array -> treat as array container return map_array_container( - partial(_apply_skew_symmetric_hybrid_diff_operator, dcoll, dq, df), - flux_matrices - ) + partial(_flux_differencing_helper, dcoll, diff_func), mats) - def _sbp_hybrid_diff_helper(diff_func, mats): - if not isinstance(mats, np.ndarray): - raise TypeError("argument must be an object array") - assert mats.dtype == object - return sum(diff_func(i, mat_i) for i, mat_i in enumerate(mats)) + assert mats.dtype == object - return _sbp_hybrid_diff_helper( - lambda i, flux_mat_i: _single_axis_hybridized_sbp_derivative_kernel( - dcoll, dq, df, i, flux_mat_i), - flux_matrices - ) + if mats.size: + sample_mat = mats[(0,)*mats.ndim] + if isinstance(sample_mat, np.ndarray): + assert sample_mat.dtype == object + # mats is an object array containing further object arrays + # -> treat as array container + return map_array_container( + partial(_flux_differencing_helper, dcoll, diff_func), mats) + + if mats.shape[-1] != dcoll.ambient_dim: + raise ValueError( + "last/innermost dimension of *mats* argument doesn't match " + "ambient dimension") + + div_result_shape = mats.shape[:-1] + + if len(div_result_shape) == 0: + return sum(diff_func(i, mat_i) for i, mat_i in enumerate(mats)) + else: + result = np.zeros(div_result_shape, dtype=object) + for idx in np.ndindex(div_result_shape): + result[idx] = sum( + diff_func(i, mat_i) for i, mat_i in enumerate(mats[idx])) + return result def volume_flux_differencing( @@ -245,9 +258,9 @@ def volume_flux_differencing( dq, df, vec: ArrayOrContainerT) -> ArrayOrContainerT: """todo. """ - def _reshape_to_numpy(shape, ary): + def _reshape(shape, ary): if not isinstance(ary, DOFArray): - return map_array_container(partial(_reshape_to_numpy, shape), ary) + return map_array_container(partial(_reshape, shape), ary) actx = ary.array_context return DOFArray( @@ -261,11 +274,11 @@ def _reshape_to_numpy(shape, ary): ) ) - return _apply_skew_symmetric_hybrid_diff_operator( + return _flux_differencing_helper( dcoll, - dq, df, - flux(_reshape_to_numpy((1, -1), vec), - _reshape_to_numpy((-1, 1), vec)) + lambda i, flux_mat_i: _single_axis_hybridized_sbp_derivative_kernel( + dcoll, dq, df, i, flux_mat_i), + flux(_reshape((1, -1), vec), _reshape((-1, 1), vec)) ) From d2f634cddc23e0d272d3e07ec8d4521d8a6e4b51 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 16 Dec 2021 23:24:01 -0600 Subject: [PATCH 14/25] Refactor flux computations --- grudge/flux_differencing.py | 21 ++------------------- grudge/models/euler.py | 27 ++++++++++++++++++++------- 2 files changed, 22 insertions(+), 26 deletions(-) diff --git a/grudge/flux_differencing.py b/grudge/flux_differencing.py index 26a2abfdf..ac0e71873 100644 --- a/grudge/flux_differencing.py +++ b/grudge/flux_differencing.py @@ -254,31 +254,14 @@ def _flux_differencing_helper(dcoll, diff_func, mats): def volume_flux_differencing( dcoll: DiscretizationCollection, - flux: Callable[[ArrayOrContainerT, ArrayOrContainerT], ArrayOrContainerT], - dq, df, vec: ArrayOrContainerT) -> ArrayOrContainerT: + dq, df, flux_matrices: ArrayOrContainerT) -> ArrayOrContainerT: """todo. """ - def _reshape(shape, ary): - if not isinstance(ary, DOFArray): - return map_array_container(partial(_reshape, shape), ary) - - actx = ary.array_context - return DOFArray( - actx, - data=tuple( - subary.reshape(grp.nelements, *shape) - for grp, subary in zip( - dcoll.discr_from_dd("vol").groups, - ary - ) - ) - ) - return _flux_differencing_helper( dcoll, lambda i, flux_mat_i: _single_axis_hybridized_sbp_derivative_kernel( dcoll, dq, df, i, flux_mat_i), - flux(_reshape((1, -1), vec), _reshape((-1, 1), vec)) + flux_matrices ) diff --git a/grudge/models/euler.py b/grudge/models/euler.py index db8a83dc2..c75ad5266 100644 --- a/grudge/models/euler.py +++ b/grudge/models/euler.py @@ -49,10 +49,14 @@ from abc import ABCMeta, abstractmethod from dataclasses import dataclass + +from functools import partial + from arraycontext import ( thaw, dataclass_array_container, - with_container_arithmetic + with_container_arithmetic, + map_array_container ) from meshmode.dof_array import DOFArray @@ -560,12 +564,21 @@ def operator(self, t, q): from functools import partial from grudge.flux_differencing import volume_flux_differencing - flux_diff = volume_flux_differencing( - dcoll, - partial(flux_chandrashekar, dcoll, gamma=gamma), - dq, df, - qtilde_allquad - ) + def _reshape(shape, ary): + if not isinstance(ary, DOFArray): + return map_array_container(partial(_reshape, shape), ary) + + return DOFArray(ary.array_context, data=tuple( + subary.reshape(grp.nelements, *shape) + # Just need group for determining the number of elements + for grp, subary in zip(dcoll.discr_from_dd("vol").groups, ary))) + + flux_matrices = flux_chandrashekar(dcoll, + _reshape((1, -1), qtilde_allquad), + _reshape((-1, 1), qtilde_allquad), + gamma=gamma) + + flux_diff = volume_flux_differencing(dcoll, dq, df, flux_matrices) def entropy_tpair(tpair): dd_intfaces = tpair.dd From 676e3a2330575992f7bcc516ac34bdeb9a63220b Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 13 Jan 2022 18:17:45 -0600 Subject: [PATCH 15/25] Update ESDG Euler operator using Euler-specific array container --- grudge/models/euler.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/grudge/models/euler.py b/grudge/models/euler.py index c75ad5266..195f89640 100644 --- a/grudge/models/euler.py +++ b/grudge/models/euler.py @@ -153,11 +153,11 @@ def conservative_to_entropy_vars(cv_state, gamma=1.4): """Converts from conserved variables (density, momentum, total energy) into entropy variables. - :arg cv_state: A :class:`EulerContainer` containing the conserved + :arg cv_state: A :class:`EulerField` containing the conserved variables. :arg gamma: The isentropic expansion factor for a single-species gas (default set to 1.4). - :returns: A :class:`EulerContainer` containing the entropy variables. + :returns: A :class:`EulerField` containing the entropy variables. """ actx = cv_state.array_context rho, u, p = conservative_to_primitive_vars(cv_state, gamma=gamma) @@ -177,11 +177,11 @@ def entropy_to_conservative_vars(ev_state, gamma=1.4): """Converts from entropy variables into conserved variables (density, momentum, total energy). - :arg ev_state: A :class:`EulerContainer` containing the entropy + :arg ev_state: A :class:`EulerField` containing the entropy variables. :arg gamma: The isentropic expansion factor for a single-species gas (default set to 1.4). - :returns: A :class:`EulerContainer` containing the conserved variables. + :returns: A :class:`EulerField` containing the conserved variables. """ actx = ev_state.array_context # See Hughes, Franca, Mallet (1986) A new finite element @@ -485,9 +485,9 @@ def ln_mean(x: DOFArray, y: DOFArray, epsilon=1e-4): + np.dot(momentum_flux, u_avg) ) - return EulerContainer(mass=mass_flux, - energy=energy_flux, - momentum=momentum_flux) + return ConservedEulerField(mass=mass_flux, + energy=energy_flux, + momentum=momentum_flux) def entropy_stable_numerical_flux_chandrashekar( From 4c9dcc98ade52ca0e7659fcc8926c898b3f10bdb Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 13 Jan 2022 23:00:51 -0600 Subject: [PATCH 16/25] Clean up modules, examples, and add convergence test --- examples/euler/sod.py | 65 ++----- examples/euler/taylor_green.py | 260 --------------------------- grudge/array_context.py | 2 - grudge/flux_differencing.py | 8 +- grudge/models/euler.py | 10 +- grudge/shortcuts.py | 118 ------------ test/test_euler_flux_differencing.py | 154 ++++++++++++++++ test/test_euler_model.py | 50 ++++++ test/test_euler_module.py | 134 -------------- 9 files changed, 221 insertions(+), 580 deletions(-) delete mode 100644 examples/euler/taylor_green.py create mode 100644 test/test_euler_flux_differencing.py delete mode 100644 test/test_euler_module.py diff --git a/examples/euler/sod.py b/examples/euler/sod.py index dda36cf71..a8699857b 100644 --- a/examples/euler/sod.py +++ b/examples/euler/sod.py @@ -1,5 +1,3 @@ -"""Minimal example of a grudge driver.""" - __copyright__ = """ Copyright (C) 2021 University of Illinois Board of Trustees """ @@ -29,17 +27,14 @@ import pyopencl.tools as cl_tools from arraycontext import thaw, freeze -from meshmode.array_context import ( - SingleGridWorkBalancingPytatoArrayContext as PytatoPyOpenCLArrayContext -) +from grudge.array_context import PytatoPyOpenCLArrayContext from grudge.models.euler import ( EntropyStableEulerOperator, - EulerContainer, + ConservedEulerField, PrescribedBC, conservative_to_primitive_vars, - primitive_to_conservative_vars ) -from grudge.shortcuts import lsrk54_step +from grudge.shortcuts import rk4_step from pytools.obj_array import make_obj_array @@ -49,31 +44,6 @@ logger = logging.getLogger(__name__) -# def sod_shock_initial_condition(nodes, t=0): -# gamma = 1.4 -# dim = len(nodes) -# x = nodes[0] -# actx = x.array_context -# zeros = 0*x - -# _x0 = 0.5 -# _rhoin = 1.0 -# _rhoout = 0.125 -# _pin = 1.0 -# _pout = 0.1 -# rhoin = zeros + _rhoin -# rhoout = zeros + _rhoout - -# x0 = zeros + _x0 -# sigma = 1e-13 -# weight = 0.5 * (1.0 - actx.np.tanh(1.0/sigma * (x - x0))) - -# rho = rhoout + (rhoin - rhoout)*weight -# p = _pout + (_pin - _pout)*weight -# u = make_obj_array([zeros for _ in range(dim)]) - -# return primitive_to_conservative_vars((rho, u, p), gamma=gamma) - def sod_shock_initial_condition(nodes, t=0): gamma = 1.4 dim = len(nodes) @@ -100,15 +70,11 @@ def sod_shock_initial_condition(nodes, t=0): energy = energyout + (energyin - energyout)*weight momentum = make_obj_array([zeros for _ in range(dim)]) - return EulerContainer(mass=mass, energy=energy, momentum=momentum) + return ConservedEulerField(mass=mass, energy=energy, momentum=momentum) -def run_sod_shock_tube(actx, - order=4, - resolution=32, - final_time=0.2, - overintegration=False, - visualize=False): +def run_sod_shock_tube( + actx, order=4, resolution=32, final_time=0.2, visualize=False): logger.info( """ @@ -116,11 +82,9 @@ def run_sod_shock_tube(actx, order: %s\n final_time: %s\n resolution: %s\n - overintegration: %s\n visualize: %s """, - order, final_time, resolution, - overintegration, visualize + order, final_time, resolution, visualize ) # eos-related parameters @@ -150,12 +114,7 @@ def run_sod_shock_tube(actx, QuadratureSimplexGroupFactory) exp_name = f"fld-sod-1d-N{order}-K{resolution}" - - if overintegration: - exp_name += "-overintegrated" - quad_tag = DISCR_TAG_QUAD - else: - quad_tag = None + quad_tag = DISCR_TAG_QUAD dcoll = DiscretizationCollection( actx, mesh, @@ -227,15 +186,14 @@ def rhs(t, q): assert norm_q < 10000 fields = thaw(freeze(fields, actx), actx) - fields = lsrk54_step(fields, t, dt, compiled_rhs) + fields = rk4_step(fields, t, dt, compiled_rhs) t += dt step += 1 # }}} -def main(ctx_factory, order=4, final_time=0.2, resolution=32, - overintegration=False, visualize=False): +def main(ctx_factory, order=4, final_time=0.2, resolution=32, visualize=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) actx = PytatoPyOpenCLArrayContext( @@ -247,7 +205,6 @@ def main(ctx_factory, order=4, final_time=0.2, resolution=32, actx, order=order, resolution=resolution, final_time=final_time, - overintegration=overintegration, visualize=visualize) @@ -258,7 +215,6 @@ def main(ctx_factory, order=4, final_time=0.2, resolution=32, parser.add_argument("--order", default=4, type=int) parser.add_argument("--tfinal", default=0.2, type=float) parser.add_argument("--resolution", default=32, type=int) - parser.add_argument("--oi", action="store_true") parser.add_argument("--visualize", action="store_true") args = parser.parse_args() @@ -267,5 +223,4 @@ def main(ctx_factory, order=4, final_time=0.2, resolution=32, order=args.order, final_time=args.tfinal, resolution=args.resolution, - overintegration=args.oi, visualize=args.visualize) diff --git a/examples/euler/taylor_green.py b/examples/euler/taylor_green.py deleted file mode 100644 index df80f644d..000000000 --- a/examples/euler/taylor_green.py +++ /dev/null @@ -1,260 +0,0 @@ -"""Minimal example of a grudge driver.""" - -__copyright__ = """ -Copyright (C) 2021 University of Illinois Board of Trustees -""" - -__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 -import pyopencl.tools as cl_tools - -from arraycontext import thaw, freeze -from grudge.array_context import PytatoPyOpenCLArrayContext -from grudge.models.euler import ( - primitive_to_conservative_vars, - conservative_to_primitive_vars, - EntropyStableEulerOperator -) -from grudge.shortcuts import lsrk54_step, lsrk144_step - -from pytools.obj_array import make_obj_array - -import grudge.op as op - -import logging -logger = logging.getLogger(__name__) - - -def tg_vortex_initial_condition(x_vec, t=0): - # Parameters chosen from https://arxiv.org/pdf/1909.12546.pdf - M = 0.05 - L = 1. - gamma = 1.4 - v0 = 1. - p0 = 1. - rho0 = gamma * (M ** 2) - - x, y, z = x_vec - actx = x.array_context - ones = 0*x + 1 - - p = p0 + rho0 * (v0 ** 2) / 16 * ( - actx.np.cos(2*x / L + actx.np.cos(2*y / L))) * actx.np.cos(2*z / L + 2) - u = v0 * actx.np.sin(x / L) * actx.np.cos(y / L) * actx.np.cos(z / L) - v = -v0 * actx.np.cos(x / L) * actx.np.sin(y / L) * actx.np.cos(z / L) - w = 0 * z - velocity = make_obj_array([u, v, w]) - rho = rho0 * ones - - return primitive_to_conservative_vars((rho, velocity, p), gamma=gamma) - - -def run_tg_vortex(actx, - order=3, - resolution=16, - cfl=1.0, - final_time=20, - dumpfreq=50, - overintegration=False, - timestepper="lsrk54", - visualize=False): - - logger.info( - """ - Taylor-Green vortex parameters:\n - order: %s\n - final_time: %s\n - resolution: %s\n - cfl: %s\n - dumpfreq: %s\n - overintegration: %s\n - timestepper: %s\n - visualize: %s - """, - order, final_time, resolution, cfl, - dumpfreq, overintegration, timestepper, visualize - ) - - # eos-related parameters - gamma = 1.4 - - # {{{ discretization - - from meshmode.mesh.generation import generate_regular_rect_mesh - - dim = 3 - box_ll = -np.pi - box_ur = np.pi - mesh = generate_regular_rect_mesh( - a=(box_ll,)*dim, - b=(box_ur,)*dim, - nelements_per_axis=(resolution,)*dim, - periodic=(True,)*dim) - - from grudge import DiscretizationCollection - from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD - from meshmode.discretization.poly_element import \ - (default_simplex_group_factory, - QuadratureSimplexGroupFactory) - - exp_name = f"fld-esdg-taylorgreen-N{order}-K{resolution}-cfl{cfl}" - - if overintegration: - exp_name += "-overintegrated" - quad_tag = DISCR_TAG_QUAD - else: - quad_tag = None - - dcoll = DiscretizationCollection( - actx, mesh, - discr_tag_to_group_factory={ - DISCR_TAG_BASE: default_simplex_group_factory(dim, order), - DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order) - } - ) - - # }}} - - # {{{ Euler operator - - euler_operator = EntropyStableEulerOperator( - dcoll, - flux_type="lf", - gamma=gamma, - quadrature_tag=quad_tag - ) - - def rhs(t, q): - return euler_operator.operator(t, q) - - compiled_rhs = actx.compile(rhs) - - fields = tg_vortex_initial_condition(thaw(dcoll.nodes(), actx)) - - if timestepper == "lsrk54": - stepper = lsrk54_step - else: - assert timestepper == "lsrk144" - stepper = lsrk144_step - - from grudge.dt_utils import h_min_from_volume - - cn = 0.5*(order + 1)**2 - dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn - - logger.info("Timestep size: %g", dt) - - # }}} - - from grudge.shortcuts import make_visualizer - - vis = make_visualizer(dcoll) - - # {{{ time stepping - - step = 0 - t = 0.0 - while t < final_time: - if step % dumpfreq == 0: - norm_q = actx.to_numpy(op.norm(dcoll, fields, 2)) - logger.info("[%04d] t = %.5f |q| = %.5e", step, t, norm_q) - if visualize: - rho, velocity, pressure = \ - conservative_to_primitive_vars(fields, gamma=gamma) - vis.write_vtk_file( - f"fld-taylor-green-vortex-{step:04d}.vtu", - [ - ("rho", rho), - ("energy", fields.energy), - ("momentum", fields.momentum), - ("velocity", velocity), - ("pressure", pressure) - ] - ) - assert norm_q < 10000 - - fields = thaw(freeze(fields, actx), actx) - fields = stepper(fields, t, dt, compiled_rhs) - t += dt - step += 1 - - # }}} - - -def main(ctx_factory, order=3, final_time=20, - resolution=16, cfl=1, - dumpfreq=50, overintegration=False, - timestepper="lsrk54", visualize=False): - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PytatoPyOpenCLArrayContext( - queue, - allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)) - ) - run_tg_vortex( - actx, - order=order, - resolution=resolution, - cfl=cfl, - final_time=final_time, - dumpfreq=dumpfreq, - overintegration=overintegration, - timestepper=timestepper, - visualize=visualize - ) - - -if __name__ == "__main__": - import argparse - - parser = argparse.ArgumentParser() - parser.add_argument("--order", default=3, type=int) - parser.add_argument("--tfinal", default=20., type=float, - help="specify final time for the simulation") - parser.add_argument("--resolution", default=8, type=int, - help="resolution in each spatial direction") - parser.add_argument("--dumpfreq", default=50, type=int, - help="step frequency for writing output or logging info") - parser.add_argument("--cfl", default=1., type=float, - help="specify cfl value to use in dt computation") - parser.add_argument("--oi", action="store_true", - help="use overintegration") - parser.add_argument("--visualize", action="store_true", - help="write out vtk output") - parser.add_argument("--timestepper", default="lsrk54", - choices=['lsrk54', 'lsrk144'], - help="specify timestepper method") - args = parser.parse_args() - - logging.basicConfig(level=logging.INFO) - main(cl.create_some_context, - order=args.order, - final_time=args.tfinal, - resolution=args.resolution, - cfl=args.cfl, - dumpfreq=args.dumpfreq, - overintegration=args.oi, - timestepper=args.timestepper, - visualize=args.visualize) diff --git a/grudge/array_context.py b/grudge/array_context.py index 2a624f863..d5d82272f 100644 --- a/grudge/array_context.py +++ b/grudge/array_context.py @@ -28,8 +28,6 @@ from meshmode.array_context import ( PyOpenCLArrayContext as _PyOpenCLArrayContextBase, - # FIXME: PytatoPyOpenCLArrayContext is literally unusable when - # using flux-differencing # TODO: Get SingleGridWorkBalancingPytatoArrayContext merged into main SingleGridWorkBalancingPytatoArrayContext as _PytatoPyOpenCLArrayContextBase, ) diff --git a/grudge/flux_differencing.py b/grudge/flux_differencing.py index ac0e71873..621b7016a 100644 --- a/grudge/flux_differencing.py +++ b/grudge/flux_differencing.py @@ -29,16 +29,14 @@ freeze ) from arraycontext.container import ArrayOrContainerT -from typing import Callable + from functools import partial + from meshmode.transform_metadata import FirstAxisIsElementsTag +from meshmode.dof_array import DOFArray from grudge.discretization import DiscretizationCollection -import grudge.dof_desc as dof_desc - -from meshmode.dof_array import DOFArray - from pytools import memoize_in, keyed_memoize_in import numpy as np diff --git a/grudge/models/euler.py b/grudge/models/euler.py index 195f89640..820cf7593 100644 --- a/grudge/models/euler.py +++ b/grudge/models/euler.py @@ -50,8 +50,6 @@ from dataclasses import dataclass -from functools import partial - from arraycontext import ( thaw, dataclass_array_container, @@ -153,11 +151,11 @@ def conservative_to_entropy_vars(cv_state, gamma=1.4): """Converts from conserved variables (density, momentum, total energy) into entropy variables. - :arg cv_state: A :class:`EulerField` containing the conserved + :arg cv_state: A :class:`ConservedEulerField` containing the conserved variables. :arg gamma: The isentropic expansion factor for a single-species gas (default set to 1.4). - :returns: A :class:`EulerField` containing the entropy variables. + :returns: A :class:`ConservedEulerField` containing the entropy variables. """ actx = cv_state.array_context rho, u, p = conservative_to_primitive_vars(cv_state, gamma=gamma) @@ -177,11 +175,11 @@ def entropy_to_conservative_vars(ev_state, gamma=1.4): """Converts from entropy variables into conserved variables (density, momentum, total energy). - :arg ev_state: A :class:`EulerField` containing the entropy + :arg ev_state: A :class:`ConservedEulerField` containing the entropy variables. :arg gamma: The isentropic expansion factor for a single-species gas (default set to 1.4). - :returns: A :class:`EulerField` containing the conserved variables. + :returns: A :class:`ConservedEulerField` containing the conserved variables. """ actx = ev_state.array_context # See Hughes, Franca, Mallet (1986) A new finite element diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py index 73f743a4e..52c77097d 100644 --- a/grudge/shortcuts.py +++ b/grudge/shortcuts.py @@ -22,124 +22,6 @@ THE SOFTWARE. """ -import numpy as np - -from dataclasses import dataclass - - -@dataclass(frozen=True) -class LSRKCoefficients: - """Dataclass which defines a given low-storage Runge-Kutta (LSRK) scheme. - The methods are determined by the provided `A`, `B` and `C` coefficient arrays. - """ - - A: np.ndarray - B: np.ndarray - C: np.ndarray - - -LSRK144NiegemannDiehlBuschCoefs = LSRKCoefficients( - A=np.array([ - 0., - -0.7188012108672410, - -0.7785331173421570, - -0.0053282796654044, - -0.8552979934029281, - -3.9564138245774565, - -1.5780575380587385, - -2.0837094552574054, - -0.7483334182761610, - -0.7032861106563359, - 0.0013917096117681, - -0.0932075369637460, - -0.9514200470875948, - -7.1151571693922548]), - B=np.array([ - 0.0367762454319673, - 0.3136296607553959, - 0.1531848691869027, - 0.0030097086818182, - 0.3326293790646110, - 0.2440251405350864, - 0.3718879239592277, - 0.6204126221582444, - 0.1524043173028741, - 0.0760894927419266, - 0.0077604214040978, - 0.0024647284755382, - 0.0780348340049386, - 5.5059777270269628]), - C=np.array([ - 0., - 0.0367762454319673, - 0.1249685262725025, - 0.2446177702277698, - 0.2476149531070420, - 0.2969311120382472, - 0.3978149645802642, - 0.5270854589440328, - 0.6981269994175695, - 0.8190890835352128, - 0.8527059887098624, - 0.8604711817462826, - 0.8627060376969976, - 0.8734213127600976])) - - -LSRK54CarpenterKennedyCoefs = LSRKCoefficients( - A=np.array([ - 0., - -567301805773/1357537059087, - -2404267990393/2016746695238, - -3550918686646/2091501179385, - -1275806237668/842570457699]), - B=np.array([ - 1432997174477/9575080441755, - 5161836677717/13612068292357, - 1720146321549/2090206949498, - 3134564353537/4481467310338, - 2277821191437/14882151754819]), - C=np.array([ - 0., - 1432997174477/9575080441755, - 2526269341429/6820363962896, - 2006345519317/3224310063776, - 2802321613138/2924317926251])) - - -EulerCoefs = LSRKCoefficients( - A=np.array([0.]), - B=np.array([1.]), - C=np.array([0.])) - - -def lsrk_step(coefs, state, t, dt, rhs): - """Take one step using a low-storage Runge-Kutta method.""" - k = 0.0 * state - for i in range(len(coefs.A)): - k = coefs.A[i]*k + dt*rhs(t + coefs.C[i]*dt, state) - state += coefs.B[i]*k - return state - - -def lsrk144_step(state, t, dt, rhs): - """Take one step using an explicit 14-stage, 4th-order, LSRK method. - - This method is derived by Niegemann, Diehl, and Busch (2012), with - an optimal stability region for advection-dominated flows. - """ - return lsrk_step(LSRK144NiegemannDiehlBuschCoefs, state, t, dt, rhs) - - -def lsrk54_step(state, t, dt, rhs): - """Take one step using an explicit 5-stage, 4th-order, LSRK method.""" - return lsrk_step(LSRK54CarpenterKennedyCoefs, state, t, dt, rhs) - - -def euler_step(state, t, dt, rhs): - """Take one step using the explicit, 1st-order accurate, Euler method.""" - return lsrk_step(EulerCoefs, state, t, dt, rhs) - def rk4_step(y, t, h, f): k1 = f(t, y) diff --git a/test/test_euler_flux_differencing.py b/test/test_euler_flux_differencing.py new file mode 100644 index 000000000..79dc6a91b --- /dev/null +++ b/test/test_euler_flux_differencing.py @@ -0,0 +1,154 @@ +__copyright__ = """ +Copyright (C) 2021 University of Illinois Board of Trustees +""" + +__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 pytest + +from grudge.array_context import PytestPytatoPyOpenCLArrayContextFactory +from arraycontext import ( + pytest_generate_tests_for_array_contexts, + thaw, freeze +) +pytest_generate_tests = pytest_generate_tests_for_array_contexts( + [PytestPytatoPyOpenCLArrayContextFactory]) + +import grudge.op as op + +import logging +logger = logging.getLogger(__name__) + + +def test_euler_esdg_vortex_convergence(actx_factory): + + from meshmode.mesh.generation import generate_regular_rect_mesh + + from grudge import DiscretizationCollection + from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD + from grudge.dt_utils import h_max_from_volume + from grudge.models.euler import ( + vortex_initial_condition, + EntropyStableEulerOperator + ) + from grudge.shortcuts import rk4_step + + from meshmode.discretization.poly_element import \ + (default_simplex_group_factory, + QuadratureSimplexGroupFactory) + + from pytools.convergence import EOCRecorder + + actx = actx_factory() + eoc_rec = EOCRecorder() + quad_tag = DISCR_TAG_QUAD + order = 2 + + for resolution in [8, 16, 32]: + + # {{{ discretization + + mesh = generate_regular_rect_mesh( + a=(0, -5), + b=(20, 5), + nelements_per_axis=(2*resolution, resolution), + periodic=(True, True)) + + discr_tag_to_group_factory = { + DISCR_TAG_BASE: default_simplex_group_factory(base_dim=2, order=order), + DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order) + } + + dcoll = DiscretizationCollection( + actx, mesh, + discr_tag_to_group_factory=discr_tag_to_group_factory + ) + h_max = actx.to_numpy(h_max_from_volume(dcoll, dim=dcoll.ambient_dim)) + nodes = thaw(dcoll.nodes(), actx) + + # }}} + + euler_operator = EntropyStableEulerOperator( + dcoll, + flux_type="lf", + gamma=1.4, + quadrature_tag=quad_tag + ) + + def rhs(t, q): + return euler_operator.operator(t, q) + + compiled_rhs = actx.compile(rhs) + + fields = vortex_initial_condition(nodes) + + from grudge.dt_utils import h_min_from_volume + + cfl = 0.125 + cn = 0.5*(order + 1)**2 + dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn + final_time = dt * 10 + + logger.info("Timestep size: %g", dt) + + # {{{ time stepping + + step = 0 + t = 0.0 + last_q = None + while t < final_time: + fields = thaw(freeze(fields, actx), actx) + fields = rk4_step(fields, t, dt, compiled_rhs) + t += dt + logger.info("[%04d] t = %.5f", step, t) + last_q = fields + last_t = t + step += 1 + + # }}} + + error_l2 = op.norm( + dcoll, + last_q - vortex_initial_condition(nodes, t=last_t), + 2 + ) + error_l2 = actx.to_numpy(error_l2) + logger.info("h_max %.5e error %.5e", h_max, error_l2) + eoc_rec.add_data_point(h_max, error_l2) + + logger.info("\n%s", eoc_rec.pretty_print(abscissa_label="h", + error_label="L2 Error")) + assert ( + eoc_rec.order_estimate() >= order + 0.5 + ) + + +# You can test individual routines by typing +# $ python test_grudge.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) + +# vim: fdm=marker diff --git a/test/test_euler_model.py b/test/test_euler_model.py index 3f5174e26..b0e5bf2a4 100644 --- a/test/test_euler_model.py +++ b/test/test_euler_model.py @@ -34,6 +34,8 @@ import grudge.op as op +import numpy as np + import logging logger = logging.getLogger(__name__) @@ -141,6 +143,54 @@ def rhs(t, q): ) +def test_entropy_variable_transformations(actx_factory): + from grudge.models.euler import ( + entropy_to_conservative_vars, + conservative_to_entropy_vars, + vortex_initial_condition + ) + + actx = actx_factory() + gamma = 1.4 # Adiabatic expansion factor for single-gas Euler model + + from meshmode.mesh.generation import generate_regular_rect_mesh + + dim = 2 + res = 5 + mesh = generate_regular_rect_mesh( + a=(0, -5), + b=(20, 5), + nelements_per_axis=(2*res, res), + periodic=(True, True)) + + from meshmode.discretization.poly_element import \ + default_simplex_group_factory + from grudge import DiscretizationCollection + from grudge.dof_desc import DISCR_TAG_BASE + + order = 3 + dcoll = DiscretizationCollection( + actx, mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: default_simplex_group_factory(dim, order) + } + ) + + # Fields in conserved variables + fields = vortex_initial_condition(thaw(dcoll.nodes(), actx)) + + # Map back and forth between entropy and conserved vars + fields_ev = conservative_to_entropy_vars(fields, gamma=gamma) + ev_fields_to_cons = entropy_to_conservative_vars(fields_ev, gamma=gamma) + + assert op.norm( + dcoll, ev_fields_to_cons.mass - fields.mass, np.inf) < 1e-13 + assert op.norm( + dcoll, ev_fields_to_cons.energy - fields.energy, np.inf) < 1e-13 + assert op.norm( + dcoll, ev_fields_to_cons.momentum - fields.momentum, np.inf) < 1e-13 + + # You can test individual routines by typing # $ python test_grudge.py 'test_routine()' diff --git a/test/test_euler_module.py b/test/test_euler_module.py deleted file mode 100644 index 79059ae9f..000000000 --- a/test/test_euler_module.py +++ /dev/null @@ -1,134 +0,0 @@ -__copyright__ = "Copyright (C) 2021 University of Illinois Board of Trustees" - -__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 - -from arraycontext import thaw - -from grudge import DiscretizationCollection -from grudge.dof_desc import DISCR_TAG_BASE - -import grudge.op as op - -from pytools.obj_array import make_obj_array - -import pytest - -from grudge.array_context import PytestPyOpenCLArrayContextFactory -from arraycontext import pytest_generate_tests_for_array_contexts -pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory]) - -import logging - -logger = logging.getLogger(__name__) - - -def test_variable_transformations(actx_factory): - from grudge.models.euler import ( - primitive_to_conservative_vars, - conservative_to_primitive_vars, - entropy_to_conservative_vars, - conservative_to_entropy_vars - ) - - actx = actx_factory() - gamma = 1.4 # Adiabatic expansion factor for single-gas Euler model - - def vortex(x_vec, t=0): - mach = 0.5 # Mach number - _x0 = 5 - epsilon = 1 # vortex strength - x, y = x_vec - actx = x.array_context - - fxyt = 1 - (((x - _x0) - t)**2 + y**2) - expterm = actx.np.exp(fxyt/2) - - u = 1 - (epsilon*y/(2*np.pi))*expterm - v = ((epsilon*(x - _x0) - t)/(2*np.pi))*expterm - - velocity = make_obj_array([u, v]) - mass = ( - 1 - ((epsilon**2*(gamma - 1)*mach**2)/(8*np.pi**2))*actx.np.exp(fxyt) - ) ** (1 / (gamma - 1)) - p = (mass**gamma)/(gamma*mach**2) - - return primitive_to_conservative_vars((mass, velocity, p), gamma=gamma) - - from meshmode.mesh.generation import generate_regular_rect_mesh - - dim = 2 - res = 5 - mesh = generate_regular_rect_mesh( - a=(0, -5), - b=(20, 5), - nelements_per_axis=(2*res, res), - periodic=(True, True)) - - from meshmode.discretization.poly_element import \ - default_simplex_group_factory - - order = 3 - dcoll = DiscretizationCollection( - actx, mesh, - discr_tag_to_group_factory={ - DISCR_TAG_BASE: default_simplex_group_factory(dim, order) - } - ) - - # Fields in conserved variables - fields = vortex(thaw(dcoll.nodes(), actx)) - - # Map back and forth between primitive and conserved vars - fields_prim = conservative_to_primitive_vars(fields, gamma=gamma) - prim_fields_to_cons = primitive_to_conservative_vars(fields_prim, gamma=gamma) - - assert op.norm( - dcoll, prim_fields_to_cons.mass - fields.mass, np.inf) < 1e-13 - assert op.norm( - dcoll, prim_fields_to_cons.energy - fields.energy, np.inf) < 1e-13 - assert op.norm( - dcoll, prim_fields_to_cons.momentum - fields.momentum, np.inf) < 1e-13 - - # Map back and forth between entropy and conserved vars - fields_ev = conservative_to_entropy_vars(fields, gamma=gamma) - ev_fields_to_cons = entropy_to_conservative_vars(fields_ev, gamma=gamma) - - assert op.norm( - dcoll, ev_fields_to_cons.mass - fields.mass, np.inf) < 1e-13 - assert op.norm( - dcoll, ev_fields_to_cons.energy - fields.energy, np.inf) < 1e-13 - assert op.norm( - dcoll, ev_fields_to_cons.momentum - fields.momentum, np.inf) < 1e-13 - - -# You can test individual routines by typing -# $ python test_grudge.py 'test_routine()' - -if __name__ == "__main__": - import sys - if len(sys.argv) > 1: - exec(sys.argv[1]) - else: - pytest.main([__file__]) From 58681b456d97707f9594abd70b5a493dbff4e669 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Tue, 18 Jan 2022 12:53:39 -0600 Subject: [PATCH 17/25] Reuse reference mass matrix routine in projection --- grudge/projection.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/grudge/projection.py b/grudge/projection.py index 0380f245f..dfbd7262a 100644 --- a/grudge/projection.py +++ b/grudge/projection.py @@ -32,8 +32,6 @@ """ -import numpy as np - from functools import partial from arraycontext import ArrayContext, map_array_container @@ -89,6 +87,7 @@ def volume_quadrature_l2_projection_matrix( vol_quad_grp.discretization_key())) def get_ref_l2_proj_mat(base_grp, vol_quad_grp): from grudge.interpolation import volume_quadrature_interpolation_matrix + from grudge.op import reference_inverse_mass_matrix vdm_q = actx.to_numpy( volume_quadrature_interpolation_matrix( @@ -96,7 +95,7 @@ def get_ref_l2_proj_mat(base_grp, vol_quad_grp): ) ) weights = vol_quad_grp.quadrature_rule().weights - inv_mass_mat = np.linalg.inv(weights * vdm_q.T @ vdm_q) + inv_mass_mat = actx.to_numpy(reference_inverse_mass_matrix(actx, base_grp)) return actx.freeze(actx.from_numpy(inv_mass_mat @ (vdm_q.T * weights))) return get_ref_l2_proj_mat(base_element_group, vol_quad_element_group) From 48cfae29aaa2c316462a8d6c45123afa7cc1f102 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Tue, 25 Jan 2022 18:57:52 -0600 Subject: [PATCH 18/25] Clean up quadrature projection routine --- grudge/flux_differencing.py | 14 +++--- grudge/op.py | 2 +- grudge/projection.py | 94 ++++++++++++++++++------------------- test/test_sbp_ops.py | 16 +++---- 4 files changed, 59 insertions(+), 67 deletions(-) diff --git a/grudge/flux_differencing.py b/grudge/flux_differencing.py index 621b7016a..402ada318 100644 --- a/grudge/flux_differencing.py +++ b/grudge/flux_differencing.py @@ -61,22 +61,20 @@ def get_skew_symetric_hybridized_diff_mats( base_grp, quad_vol_grp, face_quad_grp): from meshmode.discretization.poly_element import diff_matrices from modepy import faces_for_shape, face_normal - from grudge.projection import volume_quadrature_l2_projection_matrix from grudge.interpolation import ( volume_quadrature_interpolation_matrix, surface_quadrature_interpolation_matrix ) + from grudge.op import reference_inverse_mass_matrix # {{{ Volume operators weights = quad_vol_grp.quadrature_rule().weights vdm_q = actx.to_numpy( - volume_quadrature_interpolation_matrix(actx, base_grp, quad_vol_grp) - ) - mass_mat = weights * vdm_q.T @ vdm_q - p_mat = actx.to_numpy( - volume_quadrature_l2_projection_matrix(actx, base_grp, quad_vol_grp) - ) + volume_quadrature_interpolation_matrix(actx, base_grp, quad_vol_grp)) + inv_mass_mat = actx.to_numpy( + reference_inverse_mass_matrix(actx, base_grp)) + p_mat = inv_mass_mat @ (vdm_q.T * weights) # }}} @@ -112,7 +110,7 @@ def get_skew_symetric_hybridized_diff_mats( # {{{ Hybridized (volume + surface) operators - q_mats = [p_mat.T @ mass_mat @ diff_mat @ p_mat + q_mats = [p_mat.T @ (weights * vdm_q.T @ vdm_q) @ diff_mat @ p_mat for diff_mat in diff_matrices(base_grp)] e_mat = vf_mat @ p_mat q_skew_hybridized = np.asarray( diff --git a/grudge/op.py b/grudge/op.py index 70bd1caf7..6d1b46530 100644 --- a/grudge/op.py +++ b/grudge/op.py @@ -68,7 +68,7 @@ import grudge.dof_desc as dof_desc from grudge.interpolation import interp # noqa: F401 -from grudge.projection import project # noqa: F401 +from grudge.projection import project, volume_quadrature_project # noqa: F401 from grudge.reductions import ( # noqa: F401 norm, diff --git a/grudge/projection.py b/grudge/projection.py index dfbd7262a..e93102df5 100644 --- a/grudge/projection.py +++ b/grudge/projection.py @@ -5,6 +5,7 @@ ----------- .. autofunction:: project +.. autofunction:: volume_quadrature_project """ __copyright__ = """ @@ -31,10 +32,11 @@ THE SOFTWARE. """ +import numpy as np from functools import partial -from arraycontext import ArrayContext, map_array_container +from arraycontext import map_array_container from arraycontext.container import ArrayOrContainerT from grudge.discretization import DiscretizationCollection @@ -75,65 +77,59 @@ def project( return dcoll.connection_from_dds(src, tgt)(vec) -# {{{ Projection matrices +def volume_quadrature_project( + dcoll: DiscretizationCollection, dd_q, vec) -> ArrayOrContainerT: + """Projects a field on the quadrature discreization, described by *dd_q*, + into the polynomial space described by the volume discretization. -def volume_quadrature_l2_projection_matrix( - actx: ArrayContext, base_element_group, vol_quad_element_group): - """todo. - """ - @keyed_memoize_in( - actx, volume_quadrature_l2_projection_matrix, - lambda base_grp, vol_quad_grp: (base_grp.discretization_key(), - vol_quad_grp.discretization_key())) - def get_ref_l2_proj_mat(base_grp, vol_quad_grp): - from grudge.interpolation import volume_quadrature_interpolation_matrix - from grudge.op import reference_inverse_mass_matrix - - vdm_q = actx.to_numpy( - volume_quadrature_interpolation_matrix( - actx, base_grp, vol_quad_grp - ) - ) - weights = vol_quad_grp.quadrature_rule().weights - inv_mass_mat = actx.to_numpy(reference_inverse_mass_matrix(actx, base_grp)) - return actx.freeze(actx.from_numpy(inv_mass_mat @ (vdm_q.T * weights))) - - return get_ref_l2_proj_mat(base_element_group, vol_quad_element_group) - -# }}} - - -def volume_quadrature_project(dcoll: DiscretizationCollection, dd_q, vec): - """todo. + :arg dd_q: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. + :arg vec: a :class:`~meshmode.dof_array.DOFArray` or an + :class:`~arraycontext.container.ArrayContainer` of them. + :returns: a :class:`~meshmode.dof_array.DOFArray` or an + :class:`~arraycontext.container.ArrayContainer` like *vec*. """ if not isinstance(vec, DOFArray): return map_array_container( partial(volume_quadrature_project, dcoll, dd_q), vec ) + from grudge.geometry import area_element + from grudge.interpolation import volume_quadrature_interpolation_matrix + from grudge.op import inverse_mass + actx = vec.array_context discr = dcoll.discr_from_dd("vol") quad_discr = dcoll.discr_from_dd(dd_q) + jacobians = area_element( + actx, dcoll, dd=dd_q, + _use_geoderiv_connection=actx.supports_nonscalar_broadcasting) - return DOFArray( - actx, - data=tuple( - actx.einsum("ij,ej->ei", - volume_quadrature_l2_projection_matrix( - actx, - base_element_group=bgrp, - vol_quad_element_group=qgrp - ), - vec_i, - arg_names=("Pq_mat", "vec"), - tagged=(FirstAxisIsElementsTag(),)) - - for bgrp, qgrp, vec_i in zip( - discr.groups, - quad_discr.groups, - vec + @keyed_memoize_in( + actx, volume_quadrature_project, + lambda base_grp, vol_quad_grp: (base_grp.discretization_key(), + vol_quad_grp.discretization_key())) + def get_mat(base_grp, vol_quad_grp): + vdm_q = actx.to_numpy( + volume_quadrature_interpolation_matrix( + actx, base_grp, vol_quad_grp + ) + ) + weights = np.diag(vol_quad_grp.quadrature_rule().weights) + return actx.freeze(actx.from_numpy(vdm_q.T @ weights)) + + return inverse_mass( + dcoll, + DOFArray( + actx, + data=tuple( + actx.einsum("ij,ej,ej->ei", + get_mat(bgrp, qgrp), + jac_i, + vec_i, + arg_names=("vqw_t", "jac", "vec"), + tagged=(FirstAxisIsElementsTag(),)) + for bgrp, qgrp, vec_i, jac_i in zip( + discr.groups, quad_discr.groups, vec, jacobians) ) ) ) - -# vim: foldmethod=marker diff --git a/test/test_sbp_ops.py b/test/test_sbp_ops.py index 05e8b8ee6..793c775b2 100644 --- a/test/test_sbp_ops.py +++ b/test/test_sbp_ops.py @@ -74,11 +74,11 @@ def test_reference_element_sbp_operators(actx_factory, dim, order): from meshmode.discretization.poly_element import diff_matrices from modepy import faces_for_shape, face_normal - from grudge.projection import volume_quadrature_l2_projection_matrix from grudge.interpolation import ( volume_quadrature_interpolation_matrix, surface_quadrature_interpolation_matrix ) + from grudge.op import reference_inverse_mass_matrix for vgrp, qgrp, qfgrp in zip(volm_discr.groups, quad_discr.groups, @@ -93,18 +93,16 @@ def test_reference_element_sbp_operators(actx_factory, dim, order): weights = qgrp.quadrature_rule().weights vdm_q = actx.to_numpy( - volume_quadrature_interpolation_matrix(actx, vgrp, qgrp) - ) - mass_mat = weights * vdm_q.T @ vdm_q - p_mat = actx.to_numpy( - volume_quadrature_l2_projection_matrix(actx, vgrp, qgrp) - ) + volume_quadrature_interpolation_matrix(actx, vgrp, qgrp)) + inv_mass_mat = actx.to_numpy( + reference_inverse_mass_matrix(actx, vgrp)) + p_mat = inv_mass_mat @ (vdm_q.T * weights) # }}} # Checks Pq @ Vq = Minv @ Vq.T @ W @ Vq = I assert np.allclose(p_mat @ vdm_q, - np.identity(len(mass_mat)), rtol=1e-15) + np.identity(len(inv_mass_mat)), rtol=1e-15) # {{{ Surface operators @@ -130,7 +128,7 @@ def test_reference_element_sbp_operators(actx_factory, dim, order): # {{{ Hybridized (volume + surface) operators - q_mats = [p_mat.T @ mass_mat @ diff_mat @ p_mat + q_mats = [p_mat.T @ (weights * vdm_q.T @ vdm_q) @ diff_mat @ p_mat for diff_mat in diff_matrices(vgrp)] e_mat = vf_mat @ p_mat qtilde_mats = 0.5 * np.asarray( From d46302597799b85539738e496dca358ed9108ea7 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 27 Jan 2022 12:54:34 -0600 Subject: [PATCH 19/25] Reuse divergence helper in flux-differencing function --- grudge/flux_differencing.py | 40 ++++--------------------------------- 1 file changed, 4 insertions(+), 36 deletions(-) diff --git a/grudge/flux_differencing.py b/grudge/flux_differencing.py index 402ada318..7a73d24af 100644 --- a/grudge/flux_differencing.py +++ b/grudge/flux_differencing.py @@ -214,48 +214,16 @@ def _inv_surf_metric_deriv(): ) -def _flux_differencing_helper(dcoll, diff_func, mats): - if not isinstance(mats, np.ndarray): - # mats is not an object array -> treat as array container - return map_array_container( - partial(_flux_differencing_helper, dcoll, diff_func), mats) - - assert mats.dtype == object - - if mats.size: - sample_mat = mats[(0,)*mats.ndim] - if isinstance(sample_mat, np.ndarray): - assert sample_mat.dtype == object - # mats is an object array containing further object arrays - # -> treat as array container - return map_array_container( - partial(_flux_differencing_helper, dcoll, diff_func), mats) - - if mats.shape[-1] != dcoll.ambient_dim: - raise ValueError( - "last/innermost dimension of *mats* argument doesn't match " - "ambient dimension") - - div_result_shape = mats.shape[:-1] - - if len(div_result_shape) == 0: - return sum(diff_func(i, mat_i) for i, mat_i in enumerate(mats)) - else: - result = np.zeros(div_result_shape, dtype=object) - for idx in np.ndindex(div_result_shape): - result[idx] = sum( - diff_func(i, mat_i) for i, mat_i in enumerate(mats[idx])) - return result - - def volume_flux_differencing( dcoll: DiscretizationCollection, dq, df, flux_matrices: ArrayOrContainerT) -> ArrayOrContainerT: """todo. """ - return _flux_differencing_helper( + from grudge.op import _div_helper + + return _div_helper( dcoll, - lambda i, flux_mat_i: _single_axis_hybridized_sbp_derivative_kernel( + lambda _, i, flux_mat_i: _single_axis_hybridized_sbp_derivative_kernel( dcoll, dq, df, i, flux_mat_i), flux_matrices ) From 9ba57df378663ca65504d839789a017bf26a0d55 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 27 Jan 2022 14:01:20 -0600 Subject: [PATCH 20/25] Clean up flux-differencing module --- grudge/flux_differencing.py | 94 ++++++++++++++++++------------------- 1 file changed, 47 insertions(+), 47 deletions(-) diff --git a/grudge/flux_differencing.py b/grudge/flux_differencing.py index 7a73d24af..213df2c96 100644 --- a/grudge/flux_differencing.py +++ b/grudge/flux_differencing.py @@ -36,28 +36,25 @@ from meshmode.dof_array import DOFArray from grudge.discretization import DiscretizationCollection +from grudge.dof_desc import DOFDesc from pytools import memoize_in, keyed_memoize_in import numpy as np -def skew_symmetric_hybridized_sbp_operators( +def _reference_skew_symmetric_hybridized_sbp_operators( actx: ArrayContext, base_element_group, vol_quad_element_group, face_quad_element_group, dtype): - """todo. - """ @keyed_memoize_in( - actx, skew_symmetric_hybridized_sbp_operators, + actx, _reference_skew_symmetric_hybridized_sbp_operators, lambda base_grp, quad_vol_grp, face_quad_grp: ( base_grp.discretization_key(), quad_vol_grp.discretization_key(), - face_quad_grp.discretization_key() - ) - ) - def get_skew_symetric_hybridized_diff_mats( + face_quad_grp.discretization_key())) + def get_reference_skew_symetric_hybridized_diff_mats( base_grp, quad_vol_grp, face_quad_grp): from meshmode.discretization.poly_element import diff_matrices from modepy import faces_for_shape, face_normal @@ -100,9 +97,7 @@ def get_skew_symetric_hybridized_diff_mats( actx, base_element_group=base_grp, face_quad_element_group=face_quad_grp, - dtype=dtype - ) - ) + dtype=dtype)) zero_mat = np.zeros((nfaces*nnods_per_face, nfaces*nnods_per_face), dtype=dtype) @@ -127,19 +122,22 @@ def get_skew_symetric_hybridized_diff_mats( return actx.freeze(actx.from_numpy(q_skew_hybridized)) - return get_skew_symetric_hybridized_diff_mats( + return get_reference_skew_symetric_hybridized_diff_mats( base_element_group, vol_quad_element_group, face_quad_element_group ) -def _single_axis_hybridized_sbp_derivative_kernel( - dcoll, dq, df, xyz_axis, flux_matrix): +def _single_axis_hybridized_derivative_kernel( + dcoll, dd_quad, dd_face_quad, xyz_axis, flux_matrix): + if not dcoll._has_affine_groups(): + raise NotImplementedError("Not implemented for non-affine elements yet.") + if not isinstance(flux_matrix, DOFArray): return map_array_container( - partial(_single_axis_hybridized_sbp_derivative_kernel, - dcoll, dq, df, xyz_axis), + partial(_single_axis_hybridized_derivative_kernel, + dcoll, dd_quad, dd_face_quad, xyz_axis), flux_matrix ) @@ -150,35 +148,35 @@ def _single_axis_hybridized_sbp_derivative_kernel( volume_and_surface_quadrature_interpolation ) + actx = flux_matrix.array_context + # FIXME: This is kinda meh def inverse_jac_matrix(): - @memoize_in(dcoll, (_single_axis_hybridized_sbp_derivative_kernel, dq, df)) + @memoize_in( + dcoll, + (_single_axis_hybridized_derivative_kernel, dd_quad, dd_face_quad)) def _inv_surf_metric_deriv(): - mat = actx.np.stack( - [ - actx.np.stack( - [ - volume_and_surface_quadrature_interpolation( - dcoll, dq, df, - area_element(actx, dcoll) - * inverse_surface_metric_derivative( - actx, dcoll, - rst_ax, xyz_ax - ) - ) for rst_ax in range(dcoll.dim) - ] - ) for xyz_ax in range(dcoll.ambient_dim) - ] + return freeze( + actx.np.stack( + [ + actx.np.stack( + [ + volume_and_surface_quadrature_interpolation( + dcoll, dd_quad, dd_face_quad, + area_element(actx, dcoll) + * inverse_surface_metric_derivative( + actx, dcoll, + rst_ax, xyz_axis + ) + ) for rst_ax in range(dcoll.dim) + ] + ) for xyz_axis in range(dcoll.ambient_dim) + ] + ), + actx ) - return freeze(mat, actx) return _inv_surf_metric_deriv() - actx = flux_matrix.array_context - discr = dcoll.discr_from_dd("vol") - quad_volm_discr = dcoll.discr_from_dd(dq) - quad_face_discr = dcoll.discr_from_dd(df) - - # FIXME: This assumes affine-geometry return DOFArray( actx, data=tuple( @@ -192,7 +190,7 @@ def _inv_surf_metric_deriv(): dtype=fmat_i.dtype ), ijm_i[xyz_axis], - skew_symmetric_hybridized_sbp_operators( + _reference_skew_symmetric_hybridized_sbp_operators( actx, bgrp, qvgrp, @@ -200,13 +198,13 @@ def _inv_surf_metric_deriv(): fmat_i.dtype ), fmat_i, - arg_names=("Vh_mat_t", "inv_jac_t", "ref_Q_mat", "F_mat"), + arg_names=("Vh_mat_t", "inv_jac_t", "Q_mat", "F_mat"), tagged=(FirstAxisIsElementsTag(),)) for bgrp, qvgrp, qafgrp, fmat_i, ijm_i in zip( - discr.groups, - quad_volm_discr.groups, - quad_face_discr.groups, + dcoll.discr_from_dd("vol").groups, + dcoll.discr_from_dd(dd_quad).groups, + dcoll.discr_from_dd(dd_face_quad).groups, flux_matrix, inverse_jac_matrix() ) @@ -216,15 +214,17 @@ def _inv_surf_metric_deriv(): def volume_flux_differencing( dcoll: DiscretizationCollection, - dq, df, flux_matrices: ArrayOrContainerT) -> ArrayOrContainerT: + dd_quad: DOFDesc, + dd_face_quad: DOFDesc, + flux_matrices: ArrayOrContainerT) -> ArrayOrContainerT: """todo. """ from grudge.op import _div_helper return _div_helper( dcoll, - lambda _, i, flux_mat_i: _single_axis_hybridized_sbp_derivative_kernel( - dcoll, dq, df, i, flux_mat_i), + lambda _, i, flux_mat_i: _single_axis_hybridized_derivative_kernel( + dcoll, dd_quad, dd_face_quad, i, flux_mat_i), flux_matrices ) From f9142f4f10d8040a728d36f69de52542f8d7afee Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 27 Jan 2022 14:38:49 -0600 Subject: [PATCH 21/25] Clean up interpolation operators --- grudge/flux_differencing.py | 6 +-- grudge/interpolation.py | 85 ++++++++++--------------------------- test/test_sbp_ops.py | 4 +- 3 files changed, 26 insertions(+), 69 deletions(-) diff --git a/grudge/flux_differencing.py b/grudge/flux_differencing.py index 213df2c96..c5aa44da3 100644 --- a/grudge/flux_differencing.py +++ b/grudge/flux_differencing.py @@ -96,8 +96,7 @@ def get_reference_skew_symetric_hybridized_diff_mats( surface_quadrature_interpolation_matrix( actx, base_element_group=base_grp, - face_quad_element_group=face_quad_grp, - dtype=dtype)) + face_quad_element_group=face_quad_grp)) zero_mat = np.zeros((nfaces*nnods_per_face, nfaces*nnods_per_face), dtype=dtype) @@ -186,8 +185,7 @@ def _inv_surf_metric_deriv(): actx, base_element_group=bgrp, vol_quad_element_group=qvgrp, - face_quad_element_group=qafgrp, - dtype=fmat_i.dtype + face_quad_element_group=qafgrp ), ijm_i[xyz_axis], _reference_skew_symmetric_hybridized_sbp_operators( diff --git a/grudge/interpolation.py b/grudge/interpolation.py index 754f76abf..66a9b7ef0 100644 --- a/grudge/interpolation.py +++ b/grudge/interpolation.py @@ -36,13 +36,16 @@ from arraycontext import ( ArrayContext, - map_array_container, - thaw + map_array_container ) +from arraycontext.container import ArrayOrContainerT + from functools import partial + from meshmode.transform_metadata import FirstAxisIsElementsTag from grudge.discretization import DiscretizationCollection +from grudge.dof_desc import DOFDesc from meshmode.dof_array import DOFArray @@ -66,8 +69,6 @@ def interp(dcoll: DiscretizationCollection, src, tgt, vec): def volume_quadrature_interpolation_matrix( actx: ArrayContext, base_element_group, vol_quad_element_group): - """todo. - """ @keyed_memoize_in( actx, volume_quadrature_interpolation_matrix, lambda base_grp, vol_quad_grp: (base_grp.discretization_key(), @@ -85,10 +86,7 @@ def get_volume_vand(base_grp, vol_quad_grp): def surface_quadrature_interpolation_matrix( - actx: ArrayContext, - base_element_group, face_quad_element_group, dtype): - """todo. - """ + actx: ArrayContext, base_element_group, face_quad_element_group): @keyed_memoize_in( actx, surface_quadrature_interpolation_matrix, lambda base_grp, face_quad_grp: (base_grp.discretization_key(), @@ -121,11 +119,7 @@ def get_surface_vand(base_grp, face_quad_grp): def volume_and_surface_interpolation_matrix( actx: ArrayContext, - base_element_group, - vol_quad_element_group, - face_quad_element_group, dtype): - """todo. - """ + base_element_group, vol_quad_element_group, face_quad_element_group): @keyed_memoize_in( actx, volume_and_surface_interpolation_matrix, lambda base_grp, vol_quad_grp, face_quad_grp: ( @@ -134,21 +128,15 @@ def volume_and_surface_interpolation_matrix( face_quad_grp.discretization_key())) def get_vol_surf_interpolation_matrix(base_grp, vol_quad_grp, face_quad_grp): vq_mat = actx.to_numpy( - thaw( - volume_quadrature_interpolation_matrix( - actx, base_grp, vol_quad_grp - ), - actx - ) - ) + volume_quadrature_interpolation_matrix( + actx, + base_element_group=base_grp, + vol_quad_element_group=vol_quad_grp)) vf_mat = actx.to_numpy( - thaw( - surface_quadrature_interpolation_matrix( - actx, base_grp, face_quad_grp, dtype - ), - actx - ) - ) + surface_quadrature_interpolation_matrix( + actx, + base_element_group=base_grp, + face_quad_element_group=face_quad_grp)) return actx.freeze(actx.from_numpy(np.block([[vq_mat], [vf_mat]]))) return get_vol_surf_interpolation_matrix( @@ -158,49 +146,23 @@ def get_vol_surf_interpolation_matrix(base_grp, vol_quad_grp, face_quad_grp): # }}} -def volume_quadrature_interpolation(dcoll: DiscretizationCollection, dq, vec): - if not isinstance(vec, DOFArray): - return map_array_container( - partial(volume_quadrature_interpolation, dcoll, dq), vec - ) - - actx = vec.array_context - discr = dcoll.discr_from_dd("vol") - quad_discr = dcoll.discr_from_dd(dq) - - return DOFArray( - actx, - data=tuple( - actx.einsum("ij,ej->ei", - volume_quadrature_interpolation_matrix( - actx, - base_element_group=bgrp, - vol_quad_element_group=qgrp - ), - vec_i, - arg_names=("Vq_mat", "vec"), - tagged=(FirstAxisIsElementsTag(),)) - - for bgrp, qgrp, vec_i in zip(discr.groups, quad_discr.groups, vec) - ) - ) - - def volume_and_surface_quadrature_interpolation( - dcoll: DiscretizationCollection, dq, df, vec): + dcoll: DiscretizationCollection, + dd_quad: DOFDesc, + dd_face_quad: DOFDesc, + vec: ArrayOrContainerT) -> ArrayOrContainerT: """todo. """ if not isinstance(vec, DOFArray): return map_array_container( partial(volume_and_surface_quadrature_interpolation, - dcoll, dq, df), vec + dcoll, dd_quad, dd_face_quad), vec ) actx = vec.array_context - dtype = vec.entry_dtype discr = dcoll.discr_from_dd("vol") - quad_volm_discr = dcoll.discr_from_dd(dq) - quad_face_discr = dcoll.discr_from_dd(df) + quad_volm_discr = dcoll.discr_from_dd(dd_quad) + quad_face_discr = dcoll.discr_from_dd(dd_face_quad) return DOFArray( actx, @@ -210,8 +172,7 @@ def volume_and_surface_quadrature_interpolation( actx, base_element_group=bgrp, vol_quad_element_group=qvgrp, - face_quad_element_group=qfgrp, - dtype=dtype + face_quad_element_group=qfgrp ), vec_i, arg_names=("Vh_mat", "vec"), diff --git a/test/test_sbp_ops.py b/test/test_sbp_ops.py index 793c775b2..41bf0d6d1 100644 --- a/test/test_sbp_ops.py +++ b/test/test_sbp_ops.py @@ -70,7 +70,6 @@ def test_reference_element_sbp_operators(actx_factory, dim, order): volm_discr = dcoll.discr_from_dd("vol") quad_discr = dcoll.discr_from_dd(dd_q) quad_face_discr = dcoll.discr_from_dd(dd_f) - dtype = volm_discr.zeros(actx).entry_dtype from meshmode.discretization.poly_element import diff_matrices from modepy import faces_for_shape, face_normal @@ -119,8 +118,7 @@ def test_reference_element_sbp_operators(actx_factory, dim, order): surface_quadrature_interpolation_matrix( actx, base_element_group=vgrp, - face_quad_element_group=qfgrp, - dtype=dtype + face_quad_element_group=qfgrp ) ) From b6206afe63500550ff50294b013ba3c6a5d28d4d Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 27 Jan 2022 16:54:34 -0600 Subject: [PATCH 22/25] Clean up boundary condition application --- grudge/models/euler.py | 91 ++++++++++++++++++++++-------------------- 1 file changed, 47 insertions(+), 44 deletions(-) diff --git a/grudge/models/euler.py b/grudge/models/euler.py index 820cf7593..348de9bd2 100644 --- a/grudge/models/euler.py +++ b/grudge/models/euler.py @@ -234,7 +234,7 @@ def boundary_tpair( self, dcoll: DiscretizationCollection, dd_bc: DOFDesc, - state: ConservedEulerField, t=0): + restricted_state: ConservedEulerField, t=0): pass @@ -244,13 +244,12 @@ def boundary_tpair( self, dcoll: DiscretizationCollection, dd_bc: DOFDesc, - state: ConservedEulerField, t=0): - actx = state.array_context - dd_base = as_dofdesc("vol").with_discr_tag(DISCR_TAG_BASE) + restricted_state: ConservedEulerField, t=0): + actx = restricted_state.array_context return TracePair( dd_bc, - interior=op.project(dcoll, dd_base, dd_bc, state), + interior=restricted_state, exterior=self.prescribed_state(thaw(dcoll.nodes(dd_bc), actx), t=t) ) @@ -261,11 +260,10 @@ def boundary_tpair( self, dcoll: DiscretizationCollection, dd_bc: DOFDesc, - state: ConservedEulerField, t=0): - actx = state.array_context - dd_base = as_dofdesc("vol").with_discr_tag(DISCR_TAG_BASE) + restricted_state: ConservedEulerField, t=0): + actx = restricted_state.array_context nhat = thaw(dcoll.normal(dd_bc), actx) - interior = op.project(dcoll, dd_base, dd_bc, state) + interior = restricted_state return TracePair( dd_bc, @@ -376,13 +374,12 @@ def operator(self, t, q): dcoll = self.dcoll gamma = self.gamma qtag = self.qtag + + dd_base = DOFDesc("vol", DISCR_TAG_BASE) dq = DOFDesc("vol", qtag) df = DOFDesc("all_faces", qtag) df_int = DOFDesc("int_faces", qtag) - def interp_to_quad(u): - return op.project(dcoll, "vol", dq, u) - def interp_to_quad_surf(u): return TracePair( df_int, @@ -393,7 +390,8 @@ def interp_to_quad_surf(u): # Compute volume fluxes volume_fluxes = op.weak_local_div( dcoll, dq, - interp_to_quad(euler_volume_flux(dcoll, q, gamma=gamma)) + euler_volume_flux( + dcoll, op.project(dcoll, dd_base, dq, q), gamma=gamma) ) # Compute interior interface fluxes @@ -410,20 +408,21 @@ def interp_to_quad_surf(u): # Compute boundary fluxes if self.bdry_conditions is not None: - bc_fluxes = sum( - euler_numerical_flux( + for btag in self.bdry_conditions: + boundary_condition = self.bdry_conditions[btag] + dd_bc = as_dofdesc(btag).with_discr_tag(qtag) + bc_flux = euler_numerical_flux( dcoll, - self.bdry_conditions[btag].boundary_tpair( - dcoll, - as_dofdesc(btag).with_discr_tag(qtag), - q, + boundary_condition.boundary_tpair( + dcoll=dcoll, + dd_bc=dd_bc, + restricted_state=op.project(dcoll, dd_base, dd_bc, q), t=t ), gamma=gamma, lf_stabilization=self.lf_stabilization - ) for btag in self.bdry_conditions - ) - interface_fluxes = interface_fluxes + bc_fluxes + ) + interface_fluxes = interface_fluxes + bc_flux return op.inverse_mass( dcoll, @@ -535,13 +534,15 @@ def operator(self, t, q): gamma = self.gamma qtag = self.qtag + + dd_base = DOFDesc("vol", DISCR_TAG_BASE) dq = DOFDesc("vol", qtag) df = DOFDesc("all_faces", qtag) dcoll = self.dcoll # Interpolate cv_state to vol quad grid: u_q = V_q u - q_quad = op.project(dcoll, "vol", dq, q) + q_quad = op.project(dcoll, dd_base, dq, q) # Convert to projected entropy variables: v_q = V_h P_q v(u_q) entropy_vars = conservative_to_entropy_vars(q_quad, gamma=gamma) @@ -569,7 +570,7 @@ def _reshape(shape, ary): return DOFArray(ary.array_context, data=tuple( subary.reshape(grp.nelements, *shape) # Just need group for determining the number of elements - for grp, subary in zip(dcoll.discr_from_dd("vol").groups, ary))) + for grp, subary in zip(dcoll.discr_from_dd(dd_base).groups, ary))) flux_matrices = flux_chandrashekar(dcoll, _reshape((1, -1), qtilde_allquad), @@ -579,20 +580,15 @@ def _reshape(shape, ary): flux_diff = volume_flux_differencing(dcoll, dq, df, flux_matrices) def entropy_tpair(tpair): - dd_intfaces = tpair.dd - dd_intfaces_quad = dd_intfaces.with_discr_tag(qtag) + dd = tpair.dd + dd_quad = dd.with_discr_tag(qtag) # Interpolate entropy variables to the surface quadrature grid - vtilde_tpair = \ - op.project(dcoll, dd_intfaces, dd_intfaces_quad, tpair) + ev_tpair = op.project(dcoll, dd, dd_quad, tpair) return TracePair( - dd_intfaces_quad, + dd_quad, # Convert interior and exterior states to conserved variables - interior=entropy_to_conservative_vars( - vtilde_tpair.int, gamma=gamma - ), - exterior=entropy_to_conservative_vars( - vtilde_tpair.ext, gamma=gamma - ) + interior=entropy_to_conservative_vars(ev_tpair.int, gamma=gamma), + exterior=entropy_to_conservative_vars(ev_tpair.ext, gamma=gamma) ) # Computing interface numerical fluxes @@ -612,20 +608,27 @@ def entropy_tpair(tpair): # Compute boundary fluxes if self.bdry_conditions is not None: - bc_fluxes = sum( - entropy_stable_numerical_flux_chandrashekar( + for btag in self.bdry_conditions: + boundary_condition = self.bdry_conditions[btag] + dd_bc = as_dofdesc(btag).with_discr_tag(qtag) + bc_flux = entropy_stable_numerical_flux_chandrashekar( dcoll, - self.bdry_conditions[btag].boundary_tpair( - dcoll, - as_dofdesc(btag).with_discr_tag(qtag), - q, + boundary_condition.boundary_tpair( + dcoll=dcoll, + dd_bc=dd_bc, + # Pass modified conserved state to be used as + # the "interior" state for computing the boundary + # trace pair + restricted_state=entropy_to_conservative_vars( + op.project(dcoll, dd_base, dd_bc, proj_entropy_vars), + gamma=gamma + ), t=t ), gamma=gamma, lf_stabilization=self.lf_stabilization - ) for btag in self.bdry_conditions - ) - interface_fluxes = interface_fluxes + bc_fluxes + ) + interface_fluxes = interface_fluxes + bc_flux return op.inverse_mass( dcoll, From a67fd1e566f1cfb5555bde01955a22d7b06b7bd1 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 27 Jan 2022 19:28:14 -0600 Subject: [PATCH 23/25] Clean and combine Euler ESDG tests --- test/test_euler_flux_differencing.py | 154 --------------------------- test/test_euler_model.py | 50 +++++---- 2 files changed, 31 insertions(+), 173 deletions(-) delete mode 100644 test/test_euler_flux_differencing.py diff --git a/test/test_euler_flux_differencing.py b/test/test_euler_flux_differencing.py deleted file mode 100644 index 79dc6a91b..000000000 --- a/test/test_euler_flux_differencing.py +++ /dev/null @@ -1,154 +0,0 @@ -__copyright__ = """ -Copyright (C) 2021 University of Illinois Board of Trustees -""" - -__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 pytest - -from grudge.array_context import PytestPytatoPyOpenCLArrayContextFactory -from arraycontext import ( - pytest_generate_tests_for_array_contexts, - thaw, freeze -) -pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPytatoPyOpenCLArrayContextFactory]) - -import grudge.op as op - -import logging -logger = logging.getLogger(__name__) - - -def test_euler_esdg_vortex_convergence(actx_factory): - - from meshmode.mesh.generation import generate_regular_rect_mesh - - from grudge import DiscretizationCollection - from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD - from grudge.dt_utils import h_max_from_volume - from grudge.models.euler import ( - vortex_initial_condition, - EntropyStableEulerOperator - ) - from grudge.shortcuts import rk4_step - - from meshmode.discretization.poly_element import \ - (default_simplex_group_factory, - QuadratureSimplexGroupFactory) - - from pytools.convergence import EOCRecorder - - actx = actx_factory() - eoc_rec = EOCRecorder() - quad_tag = DISCR_TAG_QUAD - order = 2 - - for resolution in [8, 16, 32]: - - # {{{ discretization - - mesh = generate_regular_rect_mesh( - a=(0, -5), - b=(20, 5), - nelements_per_axis=(2*resolution, resolution), - periodic=(True, True)) - - discr_tag_to_group_factory = { - DISCR_TAG_BASE: default_simplex_group_factory(base_dim=2, order=order), - DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order) - } - - dcoll = DiscretizationCollection( - actx, mesh, - discr_tag_to_group_factory=discr_tag_to_group_factory - ) - h_max = actx.to_numpy(h_max_from_volume(dcoll, dim=dcoll.ambient_dim)) - nodes = thaw(dcoll.nodes(), actx) - - # }}} - - euler_operator = EntropyStableEulerOperator( - dcoll, - flux_type="lf", - gamma=1.4, - quadrature_tag=quad_tag - ) - - def rhs(t, q): - return euler_operator.operator(t, q) - - compiled_rhs = actx.compile(rhs) - - fields = vortex_initial_condition(nodes) - - from grudge.dt_utils import h_min_from_volume - - cfl = 0.125 - cn = 0.5*(order + 1)**2 - dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn - final_time = dt * 10 - - logger.info("Timestep size: %g", dt) - - # {{{ time stepping - - step = 0 - t = 0.0 - last_q = None - while t < final_time: - fields = thaw(freeze(fields, actx), actx) - fields = rk4_step(fields, t, dt, compiled_rhs) - t += dt - logger.info("[%04d] t = %.5f", step, t) - last_q = fields - last_t = t - step += 1 - - # }}} - - error_l2 = op.norm( - dcoll, - last_q - vortex_initial_condition(nodes, t=last_t), - 2 - ) - error_l2 = actx.to_numpy(error_l2) - logger.info("h_max %.5e error %.5e", h_max, error_l2) - eoc_rec.add_data_point(h_max, error_l2) - - logger.info("\n%s", eoc_rec.pretty_print(abscissa_label="h", - error_label="L2 Error")) - assert ( - eoc_rec.order_estimate() >= order + 0.5 - ) - - -# You can test individual routines by typing -# $ python test_grudge.py 'test_routine()' - -if __name__ == "__main__": - import sys - if len(sys.argv) > 1: - exec(sys.argv[1]) - else: - pytest.main([__file__]) - -# vim: fdm=marker diff --git a/test/test_euler_model.py b/test/test_euler_model.py index b0e5bf2a4..a74814988 100644 --- a/test/test_euler_model.py +++ b/test/test_euler_model.py @@ -24,13 +24,15 @@ import pytest -from grudge.array_context import PytestPyOpenCLArrayContextFactory +from grudge.array_context import \ + PytestPyOpenCLArrayContextFactory, PytestPytatoPyOpenCLArrayContextFactory from arraycontext import ( pytest_generate_tests_for_array_contexts, thaw, freeze ) pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory]) + [PytestPyOpenCLArrayContextFactory, + PytestPytatoPyOpenCLArrayContextFactory]) import grudge.op as op @@ -41,8 +43,8 @@ @pytest.mark.parametrize("order", [1, 2, 3]) -def test_euler_vortex_convergence(actx_factory, order): - +@pytest.mark.parametrize("esdg", [False, True]) +def test_euler_vortex_convergence(actx_factory, order, esdg): from meshmode.mesh.generation import generate_regular_rect_mesh from grudge import DiscretizationCollection @@ -50,7 +52,8 @@ def test_euler_vortex_convergence(actx_factory, order): from grudge.dt_utils import h_max_from_volume from grudge.models.euler import ( vortex_initial_condition, - EulerOperator + EulerOperator, + EntropyStableEulerOperator ) from grudge.shortcuts import rk4_step @@ -63,6 +66,16 @@ def test_euler_vortex_convergence(actx_factory, order): actx = actx_factory() eoc_rec = EOCRecorder() quad_tag = DISCR_TAG_QUAD + if esdg: + operator_cls = EntropyStableEulerOperator + else: + operator_cls = EulerOperator + + if esdg and not actx.supports_nonscalar_broadcasting: + pytest.xfail( + "Flux-differencing computations requires an array context " + "that supports non-scalar broadcasting" + ) for resolution in [8, 16, 32]: @@ -88,7 +101,7 @@ def test_euler_vortex_convergence(actx_factory, order): # }}} - euler_operator = EulerOperator( + euler_operator = operator_cls( dcoll, flux_type="lf", gamma=1.4, @@ -138,12 +151,10 @@ def rhs(t, q): logger.info("\n%s", eoc_rec.pretty_print(abscissa_label="h", error_label="L2 Error")) - assert ( - eoc_rec.order_estimate() >= order + 0.5 - ) + assert eoc_rec.order_estimate() >= order + 0.5 -def test_entropy_variable_transformations(actx_factory): +def test_entropy_variable_roundtrip(actx_factory): from grudge.models.euler import ( entropy_to_conservative_vars, conservative_to_entropy_vars, @@ -180,15 +191,16 @@ def test_entropy_variable_transformations(actx_factory): fields = vortex_initial_condition(thaw(dcoll.nodes(), actx)) # Map back and forth between entropy and conserved vars - fields_ev = conservative_to_entropy_vars(fields, gamma=gamma) - ev_fields_to_cons = entropy_to_conservative_vars(fields_ev, gamma=gamma) - - assert op.norm( - dcoll, ev_fields_to_cons.mass - fields.mass, np.inf) < 1e-13 - assert op.norm( - dcoll, ev_fields_to_cons.energy - fields.energy, np.inf) < 1e-13 - assert op.norm( - dcoll, ev_fields_to_cons.momentum - fields.momentum, np.inf) < 1e-13 + fields_ev = conservative_to_entropy_vars(fields, gamma) + ev_fields_to_cons = entropy_to_conservative_vars(fields_ev, gamma) + residual = ev_fields_to_cons - fields + + assert actx.to_numpy(op.norm(dcoll, residual.mass, np.inf)) < 1e-13 + assert actx.to_numpy(op.norm(dcoll, residual.energy, np.inf)) < 1e-13 + assert ( + actx.to_numpy(op.norm(dcoll, residual.momentum[i], np.inf)) < 1e-13 + for i in range(dim) + ) # You can test individual routines by typing From 7de9f810797e000577d9593cfcf22d0356b7beee Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 27 Jan 2022 20:11:54 -0600 Subject: [PATCH 24/25] Overhaul/refactor of Euler models --- grudge/models/euler.py | 351 +++++++++++++++++++++-------------------- 1 file changed, 181 insertions(+), 170 deletions(-) diff --git a/grudge/models/euler.py b/grudge/models/euler.py index 348de9bd2..5343251a5 100644 --- a/grudge/models/euler.py +++ b/grudge/models/euler.py @@ -54,7 +54,8 @@ thaw, dataclass_array_container, with_container_arithmetic, - map_array_container + map_array_container, + outer ) from meshmode.dof_array import DOFArray @@ -127,14 +128,13 @@ def vortex_initial_condition( # {{{ Variable transformation and helper routines -def conservative_to_primitive_vars(cv_state: ConservedEulerField, gamma=1.4): +def conservative_to_primitive_vars(cv_state: ConservedEulerField, gamma: float): """Converts from conserved variables (density, momentum, total energy) into primitive variables (density, velocity, pressure). :arg cv_state: A :class:`ConservedEulerField` containing the conserved variables. - :arg gamma: The isentropic expansion factor for a single-species gas - (default set to 1.4). + :arg gamma: The isentropic expansion factor. :returns: A :class:`Tuple` containing the primitive variables: (density, velocity, pressure). """ @@ -147,7 +147,7 @@ def conservative_to_primitive_vars(cv_state: ConservedEulerField, gamma=1.4): return (rho, u, p) -def conservative_to_entropy_vars(cv_state, gamma=1.4): +def conservative_to_entropy_vars(cv_state: ConservedEulerField, gamma: float): """Converts from conserved variables (density, momentum, total energy) into entropy variables. @@ -158,7 +158,7 @@ def conservative_to_entropy_vars(cv_state, gamma=1.4): :returns: A :class:`ConservedEulerField` containing the entropy variables. """ actx = cv_state.array_context - rho, u, p = conservative_to_primitive_vars(cv_state, gamma=gamma) + rho, u, p = conservative_to_primitive_vars(cv_state, gamma) u_square = sum(v ** 2 for v in u) s = actx.np.log(p) - gamma*actx.np.log(rho) @@ -171,14 +171,13 @@ def conservative_to_entropy_vars(cv_state, gamma=1.4): ) -def entropy_to_conservative_vars(ev_state, gamma=1.4): +def entropy_to_conservative_vars(ev_state: ConservedEulerField, gamma: float): """Converts from entropy variables into conserved variables (density, momentum, total energy). :arg ev_state: A :class:`ConservedEulerField` containing the entropy variables. - :arg gamma: The isentropic expansion factor for a single-species gas - (default set to 1.4). + :arg gamma: The isentropic expansion factor. :returns: A :class:`ConservedEulerField` containing the conserved variables. """ actx = ev_state.array_context @@ -205,17 +204,16 @@ def entropy_to_conservative_vars(ev_state, gamma=1.4): ) -def compute_wavespeed(cv_state: ConservedEulerField, gamma=1.4): +def compute_wavespeed(cv_state: ConservedEulerField, gamma: float): """Computes the total translational wavespeed. :arg cv_state: A :class:`ConservedEulerField` containing the conserved variables. - :arg gamma: The isentropic expansion factor for a single-species gas - (default set to 1.4). + :arg gamma: The isentropic expansion factor. :returns: A :class:`~meshmode.dof_array.DOFArray` containing local wavespeeds. """ actx = cv_state.array_context - rho, u, p = conservative_to_primitive_vars(cv_state, gamma=gamma) + rho, u, p = conservative_to_primitive_vars(cv_state, gamma) return actx.np.sqrt(np.dot(u, u)) + actx.np.sqrt(gamma * (p / rho)) @@ -283,19 +281,17 @@ def boundary_tpair( # {{{ Euler operator def euler_volume_flux( - dcoll: DiscretizationCollection, cv_state: ConservedEulerField, gamma=1.4): + dcoll: DiscretizationCollection, + cv_state: ConservedEulerField, gamma: float): """Computes the (non-linear) volume flux for the Euler operator. :arg cv_state: A :class:`ConservedEulerField` containing the conserved variables. - :arg gamma: The isentropic expansion factor for a single-species gas - (default set to 1.4). + :arg gamma: The isentropic expansion factor. :returns: A :class:`ConservedEulerField` containing the volume fluxes. """ - from arraycontext import outer - - rho, u, p = conservative_to_primitive_vars(cv_state, gamma=gamma) + rho, u, p = conservative_to_primitive_vars(cv_state, gamma) return ConservedEulerField( mass=cv_state.momentum, @@ -306,40 +302,35 @@ def euler_volume_flux( def euler_numerical_flux( dcoll: DiscretizationCollection, tpair: TracePair, - gamma=1.4, lf_stabilization=False): + gamma: float, dissipation=False): """Computes the interface numerical flux for the Euler operator. :arg tpair: A :class:`grudge.trace_pair.TracePair` containing the conserved variables on the interior and exterior sides of element facets. - :arg gamma: The isentropic expansion factor for a single-species gas - (default set to 1.4). - :arg lf_stabilization: A boolean denoting whether to apply Lax-Friedrichs + :arg gamma: The isentropic expansion factor. + :arg dissipation: A boolean denoting whether to apply Lax-Friedrichs dissipation. :returns: A :class:`ConservedEulerField` containing the interface fluxes. """ - dd_intfaces = tpair.dd - dd_allfaces = dd_intfaces.with_dtag("all_faces") q_ll = tpair.int q_rr = tpair.ext actx = q_ll.array_context flux_tpair = TracePair( tpair.dd, - interior=euler_volume_flux(dcoll, q_ll, gamma=gamma), - exterior=euler_volume_flux(dcoll, q_rr, gamma=gamma) + interior=euler_volume_flux(dcoll, q_ll, gamma), + exterior=euler_volume_flux(dcoll, q_rr, gamma) ) num_flux = flux_tpair.avg - normal = thaw(dcoll.normal(dd_intfaces), actx) - - if lf_stabilization: - from arraycontext import outer + normal = thaw(dcoll.normal(tpair.dd), actx) + if dissipation: # Compute jump penalization parameter - lam = actx.np.maximum(compute_wavespeed(q_ll, gamma=gamma), - compute_wavespeed(q_rr, gamma=gamma)) + lam = actx.np.maximum(compute_wavespeed(q_ll, gamma), + compute_wavespeed(q_rr, gamma)) num_flux -= lam*outer(tpair.diff, normal)/2 - return op.project(dcoll, dd_intfaces, dd_allfaces, num_flux @ normal) + return num_flux @ normal class EulerOperator(HyperbolicOperator): @@ -368,41 +359,46 @@ def __init__(self, dcoll: DiscretizationCollection, def max_characteristic_velocity(self, actx, **kwargs): state = kwargs["state"] - return compute_wavespeed(state, gamma=self.gamma) + return compute_wavespeed(state, self.gamma) def operator(self, t, q): dcoll = self.dcoll gamma = self.gamma qtag = self.qtag + dissipation = self.lf_stabilization dd_base = DOFDesc("vol", DISCR_TAG_BASE) - dq = DOFDesc("vol", qtag) - df = DOFDesc("all_faces", qtag) - df_int = DOFDesc("int_faces", qtag) + dd_vol_quad = DOFDesc("vol", qtag) + dd_face_quad = DOFDesc("all_faces", qtag) - def interp_to_quad_surf(u): + def interp_to_quad_surf(tpair): + dd = tpair.dd + dd_quad = dd.with_discr_tag(qtag) return TracePair( - df_int, - interior=op.project(dcoll, "int_faces", df_int, u.int), - exterior=op.project(dcoll, "int_faces", df_int, u.ext) + dd_quad, + interior=op.project(dcoll, dd, dd_quad, tpair.int), + exterior=op.project(dcoll, dd, dd_quad, tpair.ext) ) - # Compute volume fluxes - volume_fluxes = op.weak_local_div( - dcoll, dq, + interior_trace_pairs = [ + interp_to_quad_surf(tpair) + for tpair in op.interior_trace_pairs(dcoll, q) + ] + + # Compute volume derivatives + volume_derivs = op.weak_local_div( + dcoll, dd_vol_quad, euler_volume_flux( - dcoll, op.project(dcoll, dd_base, dq, q), gamma=gamma) + dcoll, op.project(dcoll, dd_base, dd_vol_quad, q), gamma) ) # Compute interior interface fluxes interface_fluxes = ( sum( - euler_numerical_flux( - dcoll, - interp_to_quad_surf(tpair), - gamma=gamma, - lf_stabilization=self.lf_stabilization - ) for tpair in op.interior_trace_pairs(dcoll, q) + op.project(dcoll, qtpair.dd, dd_face_quad, + euler_numerical_flux(dcoll, qtpair, gamma, + dissipation=dissipation)) + for qtpair in interior_trace_pairs ) ) @@ -411,22 +407,27 @@ def interp_to_quad_surf(u): for btag in self.bdry_conditions: boundary_condition = self.bdry_conditions[btag] dd_bc = as_dofdesc(btag).with_discr_tag(qtag) - bc_flux = euler_numerical_flux( + bc_flux = op.project( dcoll, - boundary_condition.boundary_tpair( - dcoll=dcoll, - dd_bc=dd_bc, - restricted_state=op.project(dcoll, dd_base, dd_bc, q), - t=t - ), - gamma=gamma, - lf_stabilization=self.lf_stabilization + dd_bc, + dd_face_quad, + euler_numerical_flux( + dcoll, + boundary_condition.boundary_tpair( + dcoll=dcoll, + dd_bc=dd_bc, + restricted_state=op.project(dcoll, dd_base, dd_bc, q), + t=t + ), + gamma, + dissipation=dissipation + ) ) interface_fluxes = interface_fluxes + bc_flux return op.inverse_mass( dcoll, - volume_fluxes - op.face_mass(dcoll, df, interface_fluxes) + volume_derivs - op.face_mass(dcoll, dd_face_quad, interface_fluxes) ) # }}} @@ -434,21 +435,23 @@ def interp_to_quad_surf(u): # {{{ Entropy stable Euler operator -def flux_chandrashekar(dcoll: DiscretizationCollection, q_ll, q_rr, gamma=1.4): +def divergence_flux_chandrashekar( + dcoll: DiscretizationCollection, + q_left: ConservedEulerField, + q_right: ConservedEulerField, gamma: float): """Two-point volume flux based on the entropy conserving and kinetic energy preserving two-point flux in: - - Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite - Volume Schemes for Compressible Euler and Navier-Stokes Equations - [DOI](https://doi.org/10.4208/cicp.170712.010313a) + Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite + Volume Schemes for Compressible Euler and Navier-Stokes Equations: + `DOI `__. - :args q_ll: an array container for the "left" state. - :args q_rr: an array container for the "right" state. - :arg gamma: The isentropic expansion factor for a single-species gas - (default set to 1.4). + :args q_left: A :class:`ConservedEulerField` containing the "left" state. + :args q_right: A :class:`ConservedEulerField` containing the "right" state. + :arg gamma: The isentropic expansion factor. """ dim = dcoll.dim - actx = q_ll.array_context + actx = q_left.array_context def ln_mean(x: DOFArray, y: DOFArray, epsilon=1e-4): f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) @@ -458,25 +461,25 @@ def ln_mean(x: DOFArray, y: DOFArray, epsilon=1e-4): (y - x) / actx.np.log(y / x) ) - rho_ll, u_ll, p_ll = conservative_to_primitive_vars(q_ll, gamma=gamma) - rho_rr, u_rr, p_rr = conservative_to_primitive_vars(q_rr, gamma=gamma) + rho_left, u_left, p_left = conservative_to_primitive_vars(q_left, gamma) + rho_right, u_right, p_right = conservative_to_primitive_vars(q_right, gamma) - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr - specific_kin_ll = 0.5 * sum(v**2 for v in u_ll) - specific_kin_rr = 0.5 * sum(v**2 for v in u_rr) + beta_left = 0.5 * rho_left / p_left + beta_right = 0.5 * rho_right / p_right + specific_kin_left = 0.5 * sum(v**2 for v in u_left) + specific_kin_right = 0.5 * sum(v**2 for v in u_right) - rho_avg = 0.5 * (rho_ll + rho_rr) - rho_mean = ln_mean(rho_ll, rho_rr) - beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - u_avg = 0.5 * (u_ll + u_rr) + rho_avg = 0.5 * (rho_left + rho_right) + rho_mean = ln_mean(rho_left, rho_right) + beta_mean = ln_mean(beta_left, beta_right) + beta_avg = 0.5 * (beta_left + beta_right) + u_avg = 0.5 * (u_left + u_right) p_mean = 0.5 * rho_avg / beta_avg - velocity_square_avg = specific_kin_ll + specific_kin_rr + velocity_square_avg = specific_kin_left + specific_kin_right mass_flux = rho_mean * u_avg - momentum_flux = np.outer(mass_flux, u_avg) + np.eye(dim) * p_mean + momentum_flux = outer(mass_flux, u_avg) + np.eye(dim) * p_mean energy_flux = ( mass_flux * 0.5 * (1/(gamma - 1)/beta_mean - velocity_square_avg) + np.dot(momentum_flux, u_avg) @@ -489,40 +492,36 @@ def ln_mean(x: DOFArray, y: DOFArray, epsilon=1e-4): def entropy_stable_numerical_flux_chandrashekar( dcoll: DiscretizationCollection, tpair: TracePair, - gamma=1.4, lf_stabilization=False): + gamma: float, dissipation=False): """Entropy stable numerical flux based on the entropy conserving and kinetic energy preserving two-point flux in: - - Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite + Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes for Compressible Euler and Navier-Stokes Equations - [DOI](https://doi.org/10.4208/cicp.170712.010313a) - """ - dd_intfaces = tpair.dd - dd_allfaces = dd_intfaces.with_dtag("all_faces") - q_ll = tpair.int - q_rr = tpair.ext - actx = q_ll.array_context + `DOI `__. - num_flux = flux_chandrashekar(dcoll, q_ll, q_rr, gamma=gamma) - normal = thaw(dcoll.normal(dd_intfaces), actx) - - if lf_stabilization: - from arraycontext import outer - - rho_ll, u_ll, p_ll = conservative_to_primitive_vars(q_ll, gamma=gamma) - rho_rr, u_rr, p_rr = conservative_to_primitive_vars(q_rr, gamma=gamma) + :arg tpair: A :class:`grudge.trace_pair.TracePair` containing the conserved + variables on the interior and exterior sides of element facets. + :arg gamma: The isentropic expansion factor. + :arg dissipation: A boolean denoting whether to apply Lax-Friedrichs + dissipation. + :returns: A :class:`ConservedEulerField` containing the interface fluxes. + """ + q_int = tpair.int + q_ext = tpair.ext + actx = q_int.array_context - def compute_wavespeed(rho, u, p): - return ( - actx.np.sqrt(np.dot(u, u)) + actx.np.sqrt(gamma * (p / rho)) - ) + num_flux = divergence_flux_chandrashekar( + dcoll, q_left=q_int, q_right=q_ext, gamma=gamma) + normal = thaw(dcoll.normal(tpair.dd), actx) + if dissipation: # Compute jump penalization parameter - lam = actx.np.maximum(compute_wavespeed(rho_ll, u_ll, p_ll), - compute_wavespeed(rho_rr, u_rr, p_rr)) + lam = actx.np.maximum(compute_wavespeed(q_int, gamma), + compute_wavespeed(q_ext, gamma)) num_flux -= lam*outer(tpair.diff, normal)/2 - return op.project(dcoll, dd_intfaces, dd_allfaces, num_flux @ normal) + return num_flux @ normal class EntropyStableEulerOperator(EulerOperator): @@ -532,34 +531,42 @@ def operator(self, t, q): from grudge.interpolation import \ volume_and_surface_quadrature_interpolation + dcoll = self.dcoll gamma = self.gamma qtag = self.qtag + dissipation = self.lf_stabilization dd_base = DOFDesc("vol", DISCR_TAG_BASE) - dq = DOFDesc("vol", qtag) - df = DOFDesc("all_faces", qtag) - - dcoll = self.dcoll - - # Interpolate cv_state to vol quad grid: u_q = V_q u - q_quad = op.project(dcoll, dd_base, dq, q) - - # Convert to projected entropy variables: v_q = V_h P_q v(u_q) - entropy_vars = conservative_to_entropy_vars(q_quad, gamma=gamma) - proj_entropy_vars = volume_quadrature_project(dcoll, dq, entropy_vars) + dd_vol_quad = DOFDesc("vol", qtag) + dd_face_quad = DOFDesc("all_faces", qtag) + + # Convert to projected entropy variables: v_q = P_q v(u_q) + proj_entropy_vars = \ + volume_quadrature_project( + dcoll, dd_vol_quad, + conservative_to_entropy_vars( + # Interpolate state to vol quad grid: u_q = V_q u + op.project(dcoll, dd_base, dd_vol_quad, q), gamma)) + + def modified_conserved_vars_tpair(tpair): + dd = tpair.dd + dd_quad = dd.with_discr_tag(qtag) + # Interpolate entropy variables to the surface quadrature grid + ev_tpair = op.project(dcoll, dd, dd_quad, tpair) + return TracePair( + dd_quad, + # Convert interior and exterior states to conserved variables + interior=entropy_to_conservative_vars(ev_tpair.int, gamma), + exterior=entropy_to_conservative_vars(ev_tpair.ext, gamma) + ) - # Compute conserved state in terms of the (interpolated) - # projected entropy variables on the quad grid - qtilde_allquad = \ - entropy_to_conservative_vars( - # Interpolate projected entropy variables to - # volume + surface quadrature grids - volume_and_surface_quadrature_interpolation( - dcoll, dq, df, proj_entropy_vars - ), - gamma=gamma) + # Compute interior trace pairs containing the modified conserved + # variables (in terms of projected entropy variables) + interior_trace_pairs = [ + modified_conserved_vars_tpair(tpair) + for tpair in op.interior_trace_pairs(dcoll, proj_entropy_vars) + ] - # Compute volume derivatives using flux differencing from functools import partial from grudge.flux_differencing import volume_flux_differencing @@ -572,37 +579,35 @@ def _reshape(shape, ary): # Just need group for determining the number of elements for grp, subary in zip(dcoll.discr_from_dd(dd_base).groups, ary))) - flux_matrices = flux_chandrashekar(dcoll, - _reshape((1, -1), qtilde_allquad), - _reshape((-1, 1), qtilde_allquad), - gamma=gamma) + # Compute the (modified) conserved state in terms of the projected + # entropy variables on both the volume and surface nodes + qtilde_vol_and_surf = \ + entropy_to_conservative_vars( + # Interpolate projected entropy variables to + # volume + surface quadrature grids + volume_and_surface_quadrature_interpolation( + dcoll, dd_vol_quad, dd_face_quad, proj_entropy_vars), gamma) - flux_diff = volume_flux_differencing(dcoll, dq, df, flux_matrices) + # FIXME: These matrices are actually symmetric. Could make use + # of that to avoid redundant computation. + flux_matrices = divergence_flux_chandrashekar( + dcoll, + _reshape((1, -1), qtilde_vol_and_surf), + _reshape((-1, 1), qtilde_vol_and_surf), + gamma + ) - def entropy_tpair(tpair): - dd = tpair.dd - dd_quad = dd.with_discr_tag(qtag) - # Interpolate entropy variables to the surface quadrature grid - ev_tpair = op.project(dcoll, dd, dd_quad, tpair) - return TracePair( - dd_quad, - # Convert interior and exterior states to conserved variables - interior=entropy_to_conservative_vars(ev_tpair.int, gamma=gamma), - exterior=entropy_to_conservative_vars(ev_tpair.ext, gamma=gamma) - ) + # Compute volume derivatives using flux differencing + volume_derivs = -volume_flux_differencing( + dcoll, dd_vol_quad, dd_face_quad, flux_matrices) # Computing interface numerical fluxes interface_fluxes = ( sum( - entropy_stable_numerical_flux_chandrashekar( - dcoll, - entropy_tpair(tpair), - gamma=gamma, - lf_stabilization=self.lf_stabilization - # NOTE: Trace pairs consist of the projected entropy variables - # which will be converted to conserved variables in - # *entropy_tpair* - ) for tpair in op.interior_trace_pairs(dcoll, proj_entropy_vars) + op.project(dcoll, qtpair.dd, dd_face_quad, + entropy_stable_numerical_flux_chandrashekar( + dcoll, qtpair, gamma, dissipation=dissipation)) + for qtpair in interior_trace_pairs ) ) @@ -611,28 +616,34 @@ def entropy_tpair(tpair): for btag in self.bdry_conditions: boundary_condition = self.bdry_conditions[btag] dd_bc = as_dofdesc(btag).with_discr_tag(qtag) - bc_flux = entropy_stable_numerical_flux_chandrashekar( + bc_flux = op.project( dcoll, - boundary_condition.boundary_tpair( - dcoll=dcoll, - dd_bc=dd_bc, - # Pass modified conserved state to be used as - # the "interior" state for computing the boundary - # trace pair - restricted_state=entropy_to_conservative_vars( - op.project(dcoll, dd_base, dd_bc, proj_entropy_vars), - gamma=gamma + dd_bc, + dd_face_quad, + entropy_stable_numerical_flux_chandrashekar( + dcoll, + boundary_condition.boundary_tpair( + dcoll=dcoll, + dd_bc=dd_bc, + # Pass modified conserved state to be used as + # the "interior" state for computing the boundary + # trace pair + restricted_state=entropy_to_conservative_vars( + op.project( + dcoll, dd_base, dd_bc, proj_entropy_vars), + gamma + ), + t=t ), - t=t - ), - gamma=gamma, - lf_stabilization=self.lf_stabilization + gamma, + dissipation=dissipation + ) ) interface_fluxes = interface_fluxes + bc_flux return op.inverse_mass( dcoll, - -flux_diff - op.face_mass(dcoll, df, interface_fluxes) + volume_derivs - op.face_mass(dcoll, dd_face_quad, interface_fluxes) ) # }}} From 26d364b468686cdf7aef694625fee87768d3f00a Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 27 Jan 2022 22:49:14 -0600 Subject: [PATCH 25/25] Document flux differencing and related functions --- doc/operators.rst | 1 + grudge/flux_differencing.py | 37 ++++++++++++++++++++++++++++++++++++- grudge/models/euler.py | 8 ++++++++ 3 files changed, 45 insertions(+), 1 deletion(-) diff --git a/doc/operators.rst b/doc/operators.rst index 550a78023..881afd945 100644 --- a/doc/operators.rst +++ b/doc/operators.rst @@ -3,6 +3,7 @@ Discontinuous Galerkin operators .. automodule:: grudge.op .. automodule:: grudge.trace_pair +.. automodule:: grudge.flux_differencing Transfering data between discretizations diff --git a/grudge/flux_differencing.py b/grudge/flux_differencing.py index c5aa44da3..546485c38 100644 --- a/grudge/flux_differencing.py +++ b/grudge/flux_differencing.py @@ -1,3 +1,11 @@ +"""Grudge module for flux-differencing in entropy-stable DG methods + +Flux-differencing +----------------- + +.. autofunction:: volume_flux_differencing +""" + __copyright__ = """ Copyright (C) 2021 University of Illinois Board of Trustees """ @@ -215,7 +223,34 @@ def volume_flux_differencing( dd_quad: DOFDesc, dd_face_quad: DOFDesc, flux_matrices: ArrayOrContainerT) -> ArrayOrContainerT: - """todo. + r"""Computes the volume contribution of the DG divergence operator using + flux-differencing: + + .. math:: + + \mathrm{VOL} = \sum_{i=1}^{d} + \begin{bmatrix} + \mathbf{V}_q \\ \mathbf{V}_f + \end{bmatrix}^T + \left( + \left( \mathbf{Q}_{i} - \mathbf{Q}^T_{i} \right) + \circ \mathbf{F}_{i} + \right)\mathbf{1} + + where :math:`\circ` denotes the + `Hadamard product `__, + :math:`\mathbf{F}_{i}` are matrices whose entries are computed + as the evaluation of an entropy-conserving two-point flux function + (e.g. :func:`grudge.models.euler.divergence_flux_chandrashekar`) + and :math:`\mathbf{Q}_{i} - \mathbf{Q}^T_{i}` are the skew-symmetric + hybridized differentiation operators defined in (15) of + `this paper `__. + + :arg flux_matrices: a :class:`~meshmode.dof_array.DOFArray` or an + :class:`~arraycontext.container.ArrayContainer` of them containing + evaluations of two-point flux. + :returns: a :class:`~meshmode.dof_array.DOFArray` or an + :class:`~arraycontext.container.ArrayContainer` of them. """ from grudge.op import _div_helper diff --git a/grudge/models/euler.py b/grudge/models/euler.py index 5343251a5..095983535 100644 --- a/grudge/models/euler.py +++ b/grudge/models/euler.py @@ -4,6 +4,7 @@ ----------------- .. autoclass:: EulerOperator +.. autoclass:: EntropyStableEulerOperator Predefined initial conditions ----------------------------- @@ -20,6 +21,9 @@ .. autofunction:: euler_volume_flux .. autofunction:: euler_numerical_flux + +.. autofunction:: divergence_flux_chandrashekar +.. autofunction:: entropy_stable_numerical_flux_chandrashekar """ __copyright__ = """ @@ -525,6 +529,10 @@ def entropy_stable_numerical_flux_chandrashekar( class EntropyStableEulerOperator(EulerOperator): + """Discretizes the Euler equations using an entropy-stable + discontinuous Galerkin discretization as outlined in (15) + of `this paper `__. + """ def operator(self, t, q): from grudge.projection import volume_quadrature_project