diff --git a/examples/vortex-mpi.py b/examples/vortex-mpi.py index 47b69d933..5e74b8d3c 100644 --- a/examples/vortex-mpi.py +++ b/examples/vortex-mpi.py @@ -32,7 +32,7 @@ from grudge.shortcuts import make_visualizer from mirgecom.discretization import create_discretization_collection -from mirgecom.euler import euler_operator +from mirgecom.euler import euler_operator, entropy_stable_euler_operator from mirgecom.simutil import ( get_sim_timestep, generate_and_distribute_mesh, @@ -68,9 +68,12 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, - use_leap=False, use_profiling=False, casename=None, lazy=False, - rst_filename=None): +def main(actx_class, use_logmgr=True, ctx_factory=cl.create_some_context, + use_leap=False, use_profiling=False, casename=None, lazy=False, + use_overintegration=False, use_esdg=False, + rst_filename=None): + + """Drive the example.""" """Drive the example.""" cl_ctx = ctx_factory() @@ -97,6 +100,20 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, from mirgecom.simutil import get_reasonable_memory_pool alloc = get_reasonable_memory_pool(cl_ctx, queue) + if lazy: + actx = actx_class(comm, queue, mpi_base_tag=12000, allocator=alloc) + else: + actx = actx_class(comm, queue, allocator=alloc, force_device_scalars=True) + + if use_esdg and not actx.supports_nonscalar_broadcasting: + raise RuntimeError( + f"{actx} is not a suitable array context for using flux-differencing. " + "The underlying array context must be capable of performing basic " + "array broadcasting operations. Use PytatoPyOpenCLArrayContext instead." + ) + from mirgecom.simutil import get_reasonable_memory_pool + alloc = get_reasonable_memory_pool(cl_ctx, queue) + if lazy: actx = actx_class(comm, queue, mpi_base_tag=12000, allocator=alloc) else: @@ -109,7 +126,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, timestepper = RK4MethodBuilder("state") else: timestepper = rk4_step - t_final = 0.01 + t_final = 1.0 current_cfl = 1.0 current_dt = .001 current_t = 0 @@ -147,13 +164,19 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, local_mesh, global_nelements = generate_and_distribute_mesh(comm, generate_mesh) local_nelements = local_mesh.nelements + from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD + from meshmode.discretization.poly_element import \ + default_simplex_group_factory, QuadratureSimplexGroupFactory order = 3 dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) vis_timer = None - + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None if logmgr: logmgr_add_cl_device_info(logmgr, queue) logmgr_add_device_memory_usage(logmgr, queue) @@ -184,7 +207,7 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, eos = IdealSingleGas() vel = np.zeros(shape=(dim,)) orig = np.zeros(shape=(dim,)) - vel[:dim] = 1.0 + vel[0] = 1.0 initializer = Vortex2D(center=orig, velocity=vel) gas_model = GasModel(eos=eos) @@ -198,7 +221,10 @@ def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): boundaries = { BTAG_ALL: PrescribedFluidBoundary(boundary_state_func=boundary_solution) } - + if use_esdg: + operator_rhs = entropy_stable_euler_operator + else: + operator_rhs = euler_operator if rst_filename: current_t = restart_data["t"] current_step = restart_data["step"] @@ -276,11 +302,11 @@ def my_health_check(pressure, component_errors): health_error = False from mirgecom.simutil import check_naninf_local, check_range_local if check_naninf_local(dcoll, "vol", pressure) \ - or check_range_local(dcoll, "vol", pressure, .2, 1.02): + or check_range_local(dcoll, "vol", pressure, .2, 1.2): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") - exittol = .1 + exittol = 1.0 if max(component_errors) > exittol: health_error = True if rank == 0: @@ -356,7 +382,7 @@ def my_post_step(step, t, dt, state): def my_rhs(t, state): fluid_state = make_fluid_state(state, gas_model) - return euler_operator(dcoll, state=fluid_state, time=t, + return operator_rhs(dcoll, state=fluid_state, time=t, boundaries=boundaries, gas_model=gas_model) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, @@ -393,6 +419,10 @@ def my_rhs(t, state): import argparse casename = "vortex" parser = argparse.ArgumentParser(description=f"MIRGE-Com Example: {casename}") + parser.add_argument("--overintegration", action="store_true", + help="use overintegration in the RHS computations"), + parser.add_argument("--esdg", action="store_true", + help="use flux-differencing/entropy stable DG for inviscid computations.") parser.add_argument("--lazy", action="store_true", help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", @@ -419,7 +449,13 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file +# main(use_logmgr=args.log, use_overintegration=args.overintegration, +# use_esdg=args.esdg, use_leap=args.leap, use_profiling=args.profiling, +# casename=casename, rst_filename=rst_filename, actx_class=actx_class) + main(actx_class, use_logmgr=args.log, use_leap=args.leap, lazy=lazy, - use_profiling=args.profiling, casename=casename, rst_filename=rst_filename) + use_profiling=args.profiling, + use_overintegration=args.overintegration, + use_esdg=args.esdg, casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker diff --git a/mirgecom/euler.py b/mirgecom/euler.py index e583367a5..c40d6fe26 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -53,25 +53,41 @@ """ import numpy as np # noqa - +from arraycontext import map_array_container from meshmode.discretization.connection import FACE_RESTR_ALL from grudge.dof_desc import ( DD_VOLUME_ALL, VolumeDomainTag, DISCR_TAG_BASE, ) - -from mirgecom.gas_model import make_operator_fluid_states +from grudge.dof_desc import DOFDesc, as_dofdesc +from mirgecom.gas_model import ( + project_fluid_state, + make_operator_fluid_states, + make_fluid_state_trace_pairs, + make_entropy_projected_fluid_state, + conservative_to_entropy_vars, + entropy_to_conservative_vars +) from mirgecom.inviscid import ( inviscid_flux, inviscid_facial_flux_rusanov, - inviscid_flux_on_element_boundary + inviscid_flux_on_element_boundary, + entropy_conserving_flux_chandrashekar, + entropy_stable_inviscid_flux_rusanov ) +from meshmode.dof_array import DOFArray +from functools import partial +import grudge.op as op from mirgecom.operators import div_operator from mirgecom.utils import normalize_boundaries - - +from grudge.projection import volume_quadrature_project +from grudge.flux_differencing import volume_flux_differencing +from grudge.trace_pair import ( + TracePair, + interior_trace_pairs +) def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, inviscid_numerical_flux_func=inviscid_facial_flux_rusanov, quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, @@ -157,7 +173,162 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, return -div_operator(dcoll, dd_vol_quad, dd_allfaces_quad, inviscid_flux_vol, inviscid_flux_bnd) - +def entropy_stable_euler_operator( + discr, state, gas_model, boundaries, time=0.0, + inviscid_numerical_flux_func=entropy_stable_inviscid_flux_rusanov, + quadrature_tag=None): + """Compute RHS of the Euler flow equations using flux-differencing. + Parameters + ---------- + state: :class:`~mirgecom.gas_model.FluidState` + Fluid state object with the conserved state, and dependent + quantities. + boundaries + Dictionary of boundary functions, one for each valid btag + time + Time + gas_model: :class:`~mirgecom.gas_model.GasModel` + Physical gas model including equation of state, transport, + and kinetic properties as required by fluid state + quadrature_tag + An optional identifier denoting a particular quadrature + discretization to use during operator evaluations. + The default value is *None*. + Returns + ------- + :class:`mirgecom.fluid.ConservedVars` + Agglomerated object array of DOF arrays representing the RHS of the Euler + flow equations. + """ + dd_base = as_dofdesc("vol") + dd_vol = DOFDesc("vol", quadrature_tag) + dd_faces = DOFDesc("all_faces", quadrature_tag) + # NOTE: For single-gas this is just a fixed scalar. + # However, for mixtures, gamma is a DOFArray. For now, + # we are re-using gamma from here and *not* recomputing + # after applying entropy projections. It is unclear at this + # time whether it's strictly necessary or if this is good enough + gamma = gas_model.eos.gamma(state.cv, state.temperature) + state_quad = project_fluid_state(discr, "vol", dd_vol, state, gas_model) + + # Compute the projected (nodal) entropy variables + entropy_vars = volume_quadrature_project( + discr, dd_vol, + # Map to entropy variables + conservative_to_entropy_vars(gamma, state_quad)) + + modified_conserved_fluid_state = \ + make_entropy_projected_fluid_state(discr, dd_vol, dd_faces, + state, entropy_vars, gamma, gas_model) + + 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(discr.discr_from_dd("vol").groups, ary))) + + flux_matrices = entropy_conserving_flux_chandrashekar( + gas_model, + _reshape((1, -1), modified_conserved_fluid_state), + _reshape((-1, 1), modified_conserved_fluid_state)) + + # Compute volume derivatives using flux differencing + inviscid_flux_vol = \ + -volume_flux_differencing(discr, dd_vol, dd_faces, flux_matrices) + + def interp_to_surf_quad(utpair): + local_dd = utpair.dd + local_dd_quad = local_dd.with_discr_tag(quadrature_tag) + return TracePair( + local_dd_quad, + interior=op.project(discr, local_dd, local_dd_quad, utpair.int), + exterior=op.project(discr, local_dd, local_dd_quad, utpair.ext) + ) + + tseed_interior_pairs = None + if state.is_mixture: + # If this is a mixture, we need to exchange the temperature field because + # mixture pressure (used in the inviscid flux calculations) depends on + # temperature and we need to seed the temperature calculation for the + # (+) part of the partition boundary with the remote temperature data. + tseed_interior_pairs = [ + # Get the interior trace pairs onto the surface quadrature + # discretization (if any) + interp_to_surf_quad(tpair) + for tpair in interior_trace_pairs(discr, state.temperature) + ] + + def interp_to_surf_modified_conservedvars(gamma, utpair): + """Takes a trace pair containing the projected entropy variables + and converts them into conserved variables on the quadrature grid. + """ + local_dd = utpair.dd + local_dd_quad = local_dd.with_discr_tag(quadrature_tag) + # Interpolate entropy variables to the surface quadrature grid + vtilde_tpair = op.project(discr, local_dd, local_dd_quad, utpair) + if isinstance(gamma, DOFArray): + gamma = op.project(discr, dd_base, local_dd_quad, gamma) + return TracePair( + local_dd_quad, + # Convert interior and exterior states to conserved variables + interior=entropy_to_conservative_vars(gamma, vtilde_tpair.int), + exterior=entropy_to_conservative_vars(gamma, vtilde_tpair.ext) + ) + + cv_interior_pairs = [ + # Compute interior trace pairs using modified conservative + # variables on the quadrature grid + # (obtaining state from projected entropy variables) + interp_to_surf_modified_conservedvars(gamma, tpair) + for tpair in interior_trace_pairs(discr, entropy_vars) + ] + + boundary_states = { + # TODO: Use modified conserved vars as the input state? + # Would need to make an "entropy-projection" variant + # of *project_fluid_state* + btag: project_fluid_state( + discr, dd_base, + # Make sure we get the state on the quadrature grid + # restricted to the tag *btag* + as_dofdesc(btag).with_discr_tag(quadrature_tag), + state, gas_model) for btag in boundaries + } + + # Interior interface state pairs consisting of modified conservative + # variables and the corresponding temperature seeds + interior_states = make_fluid_state_trace_pairs(cv_interior_pairs, + gas_model, + tseed_interior_pairs) + + # Surface contributions + inviscid_flux_bnd = ( + + # Domain boundaries + sum(boundaries[btag].inviscid_divergence_flux( + discr, + # Make sure we get the state on the quadrature grid + # restricted to the tag *btag* + as_dofdesc(btag).with_discr_tag(quadrature_tag), + gas_model, + state_minus=boundary_states[btag], + time=time, + numerical_flux_func=inviscid_numerical_flux_func) + for btag in boundaries) + + # Interior boundaries (using entropy stable numerical flux) + + sum(inviscid_facial_flux(discr, gas_model=gas_model, state_pair=state_pair, + numerical_flux_func=inviscid_numerical_flux_func) + for state_pair in interior_states) + ) + + return op.inverse_mass( + discr, + inviscid_flux_vol - op.face_mass(discr, dd_faces, inviscid_flux_bnd) + ) # By default, run unitless NAME_TO_UNITS = { "mass": "", diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py index a2ecbb329..70d5319a0 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -43,9 +43,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import numpy as np # noqa from functools import partial -from meshmode.dof_array import DOFArray # noqa from dataclasses import dataclass from arraycontext import dataclass_array_container from mirgecom.fluid import ConservedVars @@ -365,8 +363,108 @@ def project_fluid_state(dcoll, src, tgt, state, gas_model, limiter_func=None): return make_fluid_state(cv=cv_sd, gas_model=gas_model, temperature_seed=temperature_seed, limiter_func=limiter_func, limiter_dd=tgt) +def make_entropy_projected_fluid_state( + discr, dd_vol, dd_faces, state, entropy_vars, gamma, gas_model): + from grudge.interpolation import volume_and_surface_quadrature_interpolation + # Interpolate to the volume and surface (concatenated) quadrature + # discretizations: v = [v_vol, v_surf] + ev_quad = volume_and_surface_quadrature_interpolation( + discr, dd_vol, dd_faces, entropy_vars) + + temperature_seed = None + if state.is_mixture: + temperature_seed = volume_and_surface_quadrature_interpolation( + discr, dd_vol, dd_faces, state.temperature) + gamma = volume_and_surface_quadrature_interpolation( + discr, dd_vol, dd_faces, gamma) + + # Convert back to conserved varaibles and use to make the new fluid state + cv_modified = entropy_to_conservative_vars(gamma, ev_quad) + + return make_fluid_state(cv=cv_modified, + gas_model=gas_model, + temperature_seed=temperature_seed) +def conservative_to_entropy_vars(gamma, state): + """Compute the entropy variables from conserved variables. + Converts from conserved variables (density, momentum, total energy) + into entropy variables. + Parameters + ---------- + state: :class:`~mirgecom.gas_model.FluidState` + The full fluid conserved and thermal state + Returns + ------- + ConservedVars + The entropy variables + """ + from mirgecom.fluid import make_conserved + + dim = state.dim + actx = state.array_context + + rho = state.mass_density + u = state.velocity + p = state.pressure + rho_species = state.species_mass_density + + u_square = sum(v ** 2 for v in u) + s = actx.np.log(p) - gamma*actx.np.log(rho) + rho_p = rho / p + rho_species_p = rho_species / p + + return make_conserved( + dim, + mass=((gamma - s)/(gamma - 1)) - 0.5 * rho_p * u_square, + energy=-rho_p, + momentum=rho_p * u, + species_mass=((gamma - s)/(gamma - 1)) - 0.5 * rho_species_p * u_square + ) + + +def entropy_to_conservative_vars(gamma, ev: ConservedVars): + """Compute the conserved variables from entropy variables *ev*. + Converts from entropy variables into conserved variables + (density, momentum, total energy). + Parameters + ---------- + ev: ConservedVars + The entropy variables + Returns + ------- + ConservedVars + The fluid conserved variables + """ + from mirgecom.fluid import make_conserved + + dim = ev.dim + actx = ev.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 * (gamma - 1) + v1 = ev_state.mass + v234 = ev_state.momentum + v5 = ev_state.energy + v6ns = ev_state.species_mass + + v_square = sum(v**2 for v in v234) + s = gamma - v1 + v_square/(2*v5) + s_species = gamma - v6ns + v_square/(2*v5) + iota = ((gamma - 1) / (-v5)**gamma)**(inv_gamma_minus_one) + rho_iota = iota * actx.np.exp(-s * inv_gamma_minus_one) + rho_iota_species = iota * actx.np.exp(-s_species * inv_gamma_minus_one) + + return make_conserved( + dim, + mass=-rho_iota * v5, + energy=rho_iota * (1 - v_square/(2*v5)), + momentum=rho_iota * v234, + species_mass=-rho_iota_species * v5 + ) def _getattr_ish(obj, name): if obj is None: return None diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index c245ee635..cf20d9abd 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -49,7 +49,16 @@ import grudge.op as op from mirgecom.fluid import make_conserved from mirgecom.utils import normalize_boundaries +from arraycontext import thaw, outer +from mirgecom.fluid import ( + make_conserved, + ConservedVars +) + +from meshmode.dof_array import DOFArray + +from pytools.obj_array import make_obj_array def inviscid_flux(state): r"""Compute the inviscid flux vectors from fluid conserved vars *cv*. @@ -317,7 +326,108 @@ def _boundary_flux(bdtag, boundary, state_minus_quad): return inviscid_flux_bnd +def entropy_conserving_flux_chandrashekar(gas_model, state_ll, state_rr): + """Compute the entropy conservative fluxes from states *cv_ll* and *cv_rr*. + This routine implements the 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) + Returns + ------- + :class:`~mirgecom.fluid.ConservedVars` + A CV object containing the matrix-valued two-point flux vectors + for each conservation equation. + """ + dim = state_ll.dim + actx = state_ll.array_context + gamma_ll = gas_model.eos.gamma(state_ll.cv, state_ll.temperature) + gamma_rr = gas_model.eos.gamma(state_rr.cv, state_rr.temperature) + + 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) + ) + + # Primitive variables for left and right states + rho_ll = state_ll.mass_density + u_ll = state_ll.velocity + p_ll = state_ll.pressure + rho_species_ll = state_ll.species_mass_density + + rho_rr = state_rr.mass_density + u_rr = state_rr.velocity + p_rr = state_rr.pressure + rho_species_rr = state_rr.species_mass_density + + 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) + rho_species_mean = make_obj_array( + [ln_mean(y_ll_i, y_rr_i) + for y_ll_i, y_rr_i in zip(rho_species_ll, rho_species_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 = outer(mass_flux, u_avg) + np.eye(dim) * p_mean + gamma = 0.5 * (gamma_ll + gamma_rr) + energy_flux = ( + mass_flux * 0.5 * ( + 1/(gamma - 1)/beta_mean - velocity_square_avg) + + np.dot(momentum_flux, u_avg) + ) + species_mass_flux = rho_species_mean.reshape(-1, 1) * u_avg + + return ConservedVars(mass=mass_flux, + energy=energy_flux, + momentum=momentum_flux, + species_mass=species_mass_flux) + + +def entropy_stable_inviscid_flux_rusanov(state_pair, gas_model, normal, **kwargs): + r"""Return the entropy stable inviscid numerical flux. + This facial flux routine is "entropy stable" in the sense that + it computes the flux average component of the interface fluxes + using an entropy conservative two-point flux + (e.g. :func:`entropy_conserving_flux_chandrashekar`). Additional + dissipation is imposed by penalizing the "jump" of the state across + interfaces. + Parameters + ---------- + state_pair: :class:`~grudge.trace_pair.TracePair` + Trace pair of :class:`~mirgecom.gas_model.FluidState` for the face upon + which the flux calculation is to be performed + Returns + ------- + :class:`~mirgecom.fluid.ConservedVars` + A CV object containing the scalar numerical fluxes at the input faces. + """ + from mirgecom.inviscid import entropy_conserving_flux_chandrashekar + + actx = state_pair.int.array_context + flux = entropy_conserving_flux_chandrashekar(gas_model, + state_pair.int, + state_pair.ext) + + # This calculates the local maximum eigenvalue of the flux Jacobian + # for a single component gas, i.e. the element-local max wavespeed |v| + c. + lam = actx.np.maximum(state_pair.int.wavespeed, state_pair.ext.wavespeed) + dissipation = -0.5*lam*outer(state_pair.ext.cv - state_pair.int.cv, normal) + return (flux + dissipation) @ normal def get_inviscid_timestep(dcoll, state, dd=DD_VOLUME_ALL): r"""Return node-local stable timestep estimate for an inviscid fluid. diff --git a/requirements.txt b/requirements.txt index e7d5df671..77f048c2a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -19,7 +19,7 @@ git+https://github.com/pythological/kanren.git#egg=miniKanren --editable git+https://github.com/inducer/modepy.git#egg=modepy --editable git+https://github.com/inducer/arraycontext.git#egg=arraycontext --editable git+https://github.com/kaushikcfd/meshmode.git#egg=meshmode ---editable git+https://github.com/inducer/grudge.git#egg=grudge +--editable git+https://github.com/hyperSuperCube/grudge.git@esdg-trial#egg=grudge --editable git+https://github.com/kaushikcfd/pytato.git#egg=pytato --editable git+https://github.com/pyrometheus/pyrometheus.git#egg=pyrometheus --editable git+https://github.com/illinois-ceesd/logpyle.git#egg=logpyle