diff --git a/doc/misc.rst b/doc/misc.rst index a7981ab20..13f579bf8 100644 --- a/doc/misc.rst +++ b/doc/misc.rst @@ -87,4 +87,9 @@ References .. [Giles_1988] Michael Giles (1988), Non-Reflecting Boundary Conditions for the Euler \ Equations, CFDL-TR-88-1 .. [Lachaud_2014] Jean Lachaud and Nagi Mansour (2014), Porous-Material Analysis Toolbox Based - on OpenFOAM and Applications, Journal of Thermophysics and Heat Transfer 28 2 + on OpenFOAM and Applications, Journal of Thermophysics and Heat Transfer 28 2 +.. [Chandrashekar_2013] Praveen Chandrashekar, Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes \ + for Compressible Euler and Navier-Stokes Equations, Communications in Computational Physics 14, 5 \ + `(DOI) `__ +.. [Renac_2021] Florent Renac, Entropy stable, robust and high-order DGSEM for the compressible multicomponent \ + Euler equations, Journal of Computational Physics, 445 `(DOI) `__ diff --git a/examples/autoignition-mpi.py b/examples/autoignition-mpi.py index 2de268441..a8506ceee 100644 --- a/examples/autoignition-mpi.py +++ b/examples/autoignition-mpi.py @@ -38,7 +38,8 @@ from mirgecom.simutil import ( get_sim_timestep, generate_and_distribute_mesh, - write_visfile + write_visfile, + ApplicationOptionsError ) from mirgecom.io import make_init_message from mirgecom.mpi import mpi_entry_point @@ -74,10 +75,9 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, use_logmgr=True, - use_leap=False, use_overintegration=False, +def main(actx_class, use_leap=False, use_overintegration=False, casename=None, rst_filename=None, log_dependent=True, - viscous_terms_on=False): + viscous_terms_on=False, use_esdg=False): """Drive example.""" if casename is None: casename = "mirgecom" @@ -90,7 +90,7 @@ def main(actx_class, use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -584,10 +584,16 @@ def my_post_step(step, t, dt, state): return make_obj_array([cv, fluid_state.temperature]), dt - from mirgecom.inviscid import inviscid_facial_flux_rusanov as inv_num_flux_func + from mirgecom.inviscid import ( + inviscid_facial_flux_rusanov, + entropy_stable_inviscid_facial_flux_rusanov + ) from mirgecom.gas_model import make_operator_fluid_states from mirgecom.navierstokes import ns_operator + inv_num_flux_func = entropy_stable_inviscid_facial_flux_rusanov if use_esdg \ + else inviscid_facial_flux_rusanov + fluid_operator = euler_operator if viscous_terms_on: fluid_operator = ns_operator @@ -606,8 +612,9 @@ def my_rhs(t, state): fluid_rhs = fluid_operator( dcoll, state=fluid_state, gas_model=gas_model, time=t, boundaries=boundaries, operator_states_quad=fluid_operator_states, - quadrature_tag=quadrature_tag, + quadrature_tag=quadrature_tag, use_esdg=use_esdg, inviscid_numerical_flux_func=inv_num_flux_func) + chem_rhs = eos.get_species_source_terms(cv, fluid_state.temperature) tseed_rhs = fluid_state.temperature - tseed cv_rhs = fluid_rhs + chem_rhs @@ -662,8 +669,8 @@ def my_rhs(t, state): help="turns on compressible Navier-Stokes RHS") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") + parser.add_argument("--esdg", action="store_true", + help="use flux-differencing/entropy stable DG for inviscid computations.") parser.add_argument("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -673,6 +680,13 @@ def my_rhs(t, state): args = parser.parse_args() from warnings import warn warn("Automatically turning off DV logging. MIRGE-Com Issue(578)") + + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + log_dependent = False viscous_terms_on = args.navierstokes @@ -687,9 +701,9 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_leap=args.leap, - use_overintegration=args.overintegration, - casename=casename, rst_filename=rst_filename, + main(actx_class, use_leap=args.leap, + use_overintegration=args.overintegration or args.esdg, + casename=casename, rst_filename=rst_filename, use_esdg=args.esdg, log_dependent=log_dependent, viscous_terms_on=args.navierstokes) # vim: foldmethod=marker diff --git a/examples/combozzle-mpi.py b/examples/combozzle-mpi.py index 604047045..dbd599e79 100644 --- a/examples/combozzle-mpi.py +++ b/examples/combozzle-mpi.py @@ -29,8 +29,6 @@ import numpy as np from functools import partial -from meshmode.array_context import PyOpenCLArrayContext - from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.shortcuts import make_visualizer from grudge.dof_desc import BoundaryDomainTag, DISCR_TAG_QUAD @@ -41,12 +39,14 @@ from mirgecom.euler import extract_vars_for_logging, units_for_logging from mirgecom.euler import euler_operator from mirgecom.navierstokes import ns_operator + from mirgecom.simutil import ( get_sim_timestep, generate_and_distribute_mesh, write_visfile, force_evaluation, - get_box_mesh + get_box_mesh, + ApplicationOptionsError ) from mirgecom.io import make_init_message from mirgecom.mpi import mpi_entry_point @@ -146,11 +146,10 @@ def __call__(self, x_vec, *, time=0.0): @mpi_entry_point -def main(use_logmgr=True, +def main(actx_class, rst_filename=None, use_overintegration=False, casename=None, - rst_filename=None, actx_class=PyOpenCLArrayContext, log_dependent=False, input_file=None, - force_eval=True): + force_eval=True, use_esdg=False): """Drive example.""" if casename is None: casename = "mirgecom" @@ -171,7 +170,7 @@ def main(use_logmgr=True, # {{{ Some discretization parameters - dim = 3 + dim = 2 order = 3 # - scales the size of the domain @@ -675,7 +674,7 @@ def vol_max(x): casename = f"{casename}-d{dim}p{order}e{global_nelements}n{nparts}" - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) if logmgr: @@ -1077,7 +1076,12 @@ def my_post_step(step, t, dt, state): return state, dt - from mirgecom.inviscid import inviscid_facial_flux_rusanov + from mirgecom.inviscid import ( + inviscid_facial_flux_rusanov, + entropy_stable_inviscid_facial_flux_rusanov + ) + inv_num_flux_func = entropy_stable_inviscid_facial_flux_rusanov if use_esdg \ + else inviscid_facial_flux_rusanov def dummy_pre_step(step, t, dt, state): if logmgr: @@ -1106,18 +1110,18 @@ def cfd_rhs(t, state): from mirgecom.gas_model import make_fluid_state fluid_state = make_fluid_state(cv=cv, gas_model=gas_model, temperature_seed=tseed) - fluid_operator_states = make_operator_fluid_states(dcoll, fluid_state, - gas_model, boundaries, - quadrature_tag) + fluid_operator_states = make_operator_fluid_states( + dcoll, fluid_state, gas_model, boundaries, quadrature_tag=quadrature_tag) if inviscid_only: fluid_rhs = \ euler_operator( dcoll, state=fluid_state, time=t, boundaries=boundaries, gas_model=gas_model, - inviscid_numerical_flux_func=inviscid_facial_flux_rusanov, + inviscid_numerical_flux_func=inv_num_flux_func, quadrature_tag=quadrature_tag, - operator_states_quad=fluid_operator_states) + operator_states_quad=fluid_operator_states, + use_esdg=use_esdg) else: grad_cv = grad_cv_operator(dcoll, gas_model, boundaries, fluid_state, time=t, @@ -1128,33 +1132,29 @@ def cfd_rhs(t, state): ns_operator( dcoll, state=fluid_state, time=t, boundaries=boundaries, gas_model=gas_model, quadrature_tag=quadrature_tag, - inviscid_numerical_flux_func=inviscid_facial_flux_rusanov) + inviscid_numerical_flux_func=inv_num_flux_func, + operator_states_quad=fluid_operator_states, + use_esdg=use_esdg) - if inert_only: - chem_rhs = 0*fluid_rhs - else: - chem_rhs = eos.get_species_source_terms(cv, fluid_state.temperature) + if not inert_only: + fluid_rhs = fluid_rhs + eos.get_species_source_terms( + cv, fluid_state.temperature) if av_on: alpha_f = compute_av_alpha_field(fluid_state) indicator = smoothness_indicator(dcoll, fluid_state.mass_density, kappa=kappa_sc, s0=s0_sc) - av_rhs = av_laplacian_operator( + fluid_rhs = fluid_rhs + av_laplacian_operator( dcoll, fluid_state=fluid_state, boundaries=boundaries, time=t, gas_model=gas_model, grad_cv=grad_cv, operator_states_quad=fluid_operator_states, alpha=alpha_f, s0=s0_sc, kappa=kappa_sc, indicator=indicator) - else: - av_rhs = 0*fluid_rhs if sponge_on: - sponge_rhs = _sponge(fluid_state.cv) - else: - sponge_rhs = 0*fluid_rhs + fluid_rhs = fluid_rhs + _sponge(fluid_state.cv) - fluid_rhs = fluid_rhs + chem_rhs + av_rhs + sponge_rhs - tseed_rhs = fluid_state.temperature - tseed + tseed_rhs = actx.zeros_like(fluid_state.temperature) return make_obj_array([fluid_rhs, tseed_rhs]) @@ -1245,8 +1245,8 @@ def dummy_rhs(t, state): help="Turn off force lazy eval between timesteps") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") + parser.add_argument("--esdg", action="store_true", + help="use entropy-stable for inviscid terms") parser.add_argument("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -1254,8 +1254,16 @@ def dummy_rhs(t, state): parser.add_argument("--restart_file", help="root name of restart file") parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn warn("Automatically turning off DV logging. MIRGE-Com Issue(578)") + + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + log_dependent = False force_eval = not args.no_force @@ -1279,9 +1287,9 @@ def dummy_rhs(t, state): print(f"Calling main: {time.ctime(time.time())}") - main(use_logmgr=args.log, input_file=input_file, - use_overintegration=args.overintegration, - casename=casename, rst_filename=rst_filename, actx_class=actx_class, - log_dependent=log_dependent, force_eval=force_eval) + main(actx_class, input_file=input_file, + use_overintegration=args.overintegration or args.esdg, + casename=casename, rst_filename=rst_filename, + log_dependent=log_dependent, force_eval=force_eval, use_esdg=args.esdg) # vim: foldmethod=marker diff --git a/examples/doublemach-mpi.py b/examples/doublemach-mpi.py index 0cdafdf81..6e396a577 100644 --- a/examples/doublemach-mpi.py +++ b/examples/doublemach-mpi.py @@ -47,7 +47,7 @@ from mirgecom.initializers import DoubleMachReflection from mirgecom.eos import IdealSingleGas from mirgecom.transport import SimpleTransport -from mirgecom.simutil import get_sim_timestep +from mirgecom.simutil import get_sim_timestep, ApplicationOptionsError from logpyle import set_dt from mirgecom.euler import extract_vars_for_logging, units_for_logging from mirgecom.logging_quantities import ( @@ -115,13 +115,10 @@ def get_doublemach_mesh(): @mpi_entry_point -def main(use_logmgr=True, +def main(actx_class, use_esdg=False, use_leap=False, use_overintegration=False, - casename=None, rst_filename=None, actx_class=None): + casename=None, rst_filename=None): """Drive the example.""" - if actx_class is None: - raise RuntimeError("Array context class missing.") - if casename is None: casename = "mirgecom" @@ -130,7 +127,7 @@ def main(use_logmgr=True, rank = comm.Get_rank() nparts = comm.Get_size() - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -391,7 +388,8 @@ def my_rhs(t, state): return ( euler_operator(dcoll, state=fluid_state, time=t, boundaries=boundaries, - gas_model=gas_model, quadrature_tag=quadrature_tag) + gas_model=gas_model, quadrature_tag=quadrature_tag, + use_esdg=use_esdg) + av_laplacian_operator(dcoll, fluid_state=fluid_state, boundaries=boundaries, time=t, gas_model=gas_model, @@ -437,8 +435,8 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") + parser.add_argument("--esdg", action="store_true", + help="use flux-differencing/entropy stable DG for inviscid computations.") parser.add_argument("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -447,6 +445,13 @@ def my_rhs(t, state): parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class(lazy=args.lazy, distributed=True, profiling=args.profiling, @@ -459,8 +464,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(use_logmgr=args.log, use_leap=args.leap, - use_overintegration=args.overintegration, - casename=casename, rst_filename=rst_filename, actx_class=actx_class) + main(actx_class, use_leap=args.leap, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, + casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker diff --git a/examples/doublemach_physical_av-mpi.py b/examples/doublemach_physical_av-mpi.py index f82299471..025c698ed 100644 --- a/examples/doublemach_physical_av-mpi.py +++ b/examples/doublemach_physical_av-mpi.py @@ -119,13 +119,10 @@ def get_doublemach_mesh(): @mpi_entry_point -def main(use_logmgr=True, +def main(actx_class, use_esdg=False, use_leap=False, use_overintegration=False, - casename=None, rst_filename=None, actx_class=None): + casename=None, rst_filename=None): """Drive the example.""" - if actx_class is None: - raise RuntimeError("Array context class missing.") - if casename is None: casename = "mirgecom" @@ -145,7 +142,7 @@ def main(use_logmgr=True, os.makedirs(log_dir) comm.Barrier() - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=logname, mode="wo", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -595,7 +592,8 @@ def _my_rhs(t, state): return ( euler_operator(dcoll, state=fluid_state, time=t, boundaries=boundaries, - gas_model=gas_model, quadrature_tag=quadrature_tag) + gas_model=gas_model, quadrature_tag=quadrature_tag, + use_esdg=use_esdg) ) def _my_rhs_av(t, state): @@ -605,7 +603,8 @@ def _my_rhs_av(t, state): return ( euler_operator(dcoll, state=fluid_state, time=t, boundaries=boundaries, - gas_model=gas_model, quadrature_tag=quadrature_tag) + gas_model=gas_model, quadrature_tag=quadrature_tag, + use_esdg=use_esdg) + av_laplacian_operator(dcoll, fluid_state=fluid_state, boundaries=boundaries, time=t, gas_model=gas_model, @@ -623,7 +622,8 @@ def _my_rhs_phys_visc_av(t, state): return ( ns_operator(dcoll, state=fluid_state, time=t, boundaries=boundaries, - gas_model=gas_model, quadrature_tag=quadrature_tag) + gas_model=gas_model, quadrature_tag=quadrature_tag, + use_esdg=use_esdg) ) def _my_rhs_phys_visc_div_av(t, state): @@ -647,7 +647,7 @@ def _my_rhs_phys_visc_div_av(t, state): ns_operator(dcoll, state=fluid_state, time=t, boundaries=boundaries, gas_model=gas_model, quadrature_tag=quadrature_tag, - grad_cv=grad_cv) + grad_cv=grad_cv, use_esdg=use_esdg) ) my_rhs = (_my_rhs if use_av == 0 else _my_rhs_av if use_av == 1 else @@ -708,8 +708,8 @@ def _my_rhs_phys_visc_div_av(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") + parser.add_argument("--esdg", action="store_true", + help="use flux-differencing/entropy stable DG for inviscid computations.") parser.add_argument("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -718,6 +718,14 @@ def _my_rhs_phys_visc_div_av(t, state): parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class(lazy=args.lazy, distributed=True, @@ -731,8 +739,8 @@ def _my_rhs_phys_visc_div_av(t, state): if args.restart_file: rst_filename = args.restart_file - main(use_logmgr=args.log, use_leap=args.leap, - use_overintegration=args.overintegration, - casename=casename, rst_filename=rst_filename, actx_class=actx_class) + main(actx_class, use_leap=args.leap, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, + casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker diff --git a/examples/heat-source-mpi.py b/examples/heat-source-mpi.py index 53beb954e..d2ec592fd 100644 --- a/examples/heat-source-mpi.py +++ b/examples/heat-source-mpi.py @@ -47,14 +47,15 @@ @mpi_entry_point -def main(actx_class, use_logmgr=True, +def main(actx_class, use_esdg=False, + use_overintegration=False, use_leap=False, casename=None, rst_filename=None): """Run the example.""" from mpi4py import MPI comm = MPI.COMM_WORLD num_parts = comm.Get_size() - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename="heat-source.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -99,6 +100,12 @@ def main(actx_class, use_logmgr=True, dcoll = create_discretization_collection(actx, local_mesh, order=order) + from grudge.dof_desc import DISCR_TAG_QUAD + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None # noqa + if dim == 2: # no deep meaning here, just a fudge factor dt = 0.0025/(nel_1d*order**2) @@ -137,7 +144,8 @@ def main(actx_class, use_logmgr=True, def rhs(t, u): return ( - diffusion_operator(dcoll, kappa=1, boundaries=boundaries, u=u) + diffusion_operator(dcoll, kappa=1, boundaries=boundaries, u=u, + quadrature_tag=quadrature_tag) + actx.np.exp(-np.dot(nodes, nodes)/source_width**2)) compiled_rhs = actx.compile(rhs) @@ -182,8 +190,8 @@ def rhs(t, u): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") + parser.add_argument("--overintegration", action="store_true", + help="turn on overintegration.") parser.add_argument("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -203,7 +211,8 @@ def rhs(t, u): if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_leap=args.leap, + main(actx_class, use_leap=args.leap, + use_overintegration=args.overintegration, casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker diff --git a/examples/hotplate-mpi.py b/examples/hotplate-mpi.py index a11adb9ed..6cc8e7ab2 100644 --- a/examples/hotplate-mpi.py +++ b/examples/hotplate-mpi.py @@ -72,13 +72,9 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(use_logmgr=True, - use_leap=False, casename=None, - rst_filename=None, actx_class=None): +def main(actx_class, use_esdg=False, use_leap=False, casename=None, + use_overintegration=False, rst_filename=None): """Drive the example.""" - if actx_class is None: - raise RuntimeError("Array context class missing.") - if casename is None: casename = "mirgecom" @@ -90,7 +86,7 @@ def main(use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -151,6 +147,12 @@ def main(use_logmgr=True, dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) + from grudge.dof_desc import DISCR_TAG_QUAD + 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) @@ -390,7 +392,8 @@ def my_post_step(step, t, dt, state): def my_rhs(t, state): fluid_state = make_fluid_state(state, gas_model) return ns_operator(dcoll, boundaries=boundaries, state=fluid_state, - time=t, gas_model=gas_model) + time=t, gas_model=gas_model, use_esdg=use_esdg, + quadrature_tag=quadrature_tag) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, t_final, constant_cfl) @@ -435,6 +438,10 @@ def my_rhs(t, state): help="turn on detailed performance profiling") parser.add_argument("--log", action="store_true", default=True, help="turn on logging") + 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("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -443,6 +450,14 @@ def my_rhs(t, state): parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -454,7 +469,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(use_logmgr=args.log, use_leap=args.leap, - casename=casename, rst_filename=rst_filename, actx_class=actx_class) + main(actx_class, use_leap=args.leap, + casename=casename, rst_filename=rst_filename, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg) # vim: foldmethod=marker diff --git a/examples/lump-mpi.py b/examples/lump-mpi.py index 63c7ee68c..7ea67efb7 100644 --- a/examples/lump-mpi.py +++ b/examples/lump-mpi.py @@ -27,6 +27,7 @@ import numpy as np from functools import partial +from grudge.dof_desc import DISCR_TAG_QUAD from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.shortcuts import make_visualizer @@ -42,7 +43,7 @@ from mirgecom.integrators import rk4_step from mirgecom.steppers import advance_state from mirgecom.boundary import PrescribedFluidBoundary -from mirgecom.initializers import Lump +from mirgecom.initializers import Lump, Uniform # noqa from mirgecom.eos import IdealSingleGas from logpyle import IntervalTimer, set_dt @@ -65,7 +66,8 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, use_logmgr=True, +def main(actx_class, use_esdg=False, + use_overintegration=False, use_leap=False, casename=None, rst_filename=None): """Drive example.""" if casename is None: @@ -79,7 +81,7 @@ def main(actx_class, use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -93,7 +95,7 @@ def main(actx_class, use_logmgr=True, timestepper = RK4MethodBuilder("state") else: timestepper = rk4_step - t_final = 0.01 + t_final = .005 current_cfl = 1.0 current_dt = .001 current_t = 0 @@ -136,6 +138,11 @@ def main(actx_class, use_logmgr=True, dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None + vis_timer = None if logmgr: @@ -162,6 +169,8 @@ def main(actx_class, use_logmgr=True, orig = np.zeros(shape=(dim,)) vel[:dim] = 1.0 initializer = Lump(dim=dim, center=orig, velocity=vel) + # initializer = Uniform(dim=dim, velocity=vel) + from mirgecom.gas_model import GasModel, make_fluid_state gas_model = GasModel(eos=eos) @@ -235,11 +244,14 @@ def my_write_restart(state, step, t): def my_health_check(dv, state, exact): health_error = False from mirgecom.simutil import check_naninf_local, check_range_local - if check_naninf_local(dcoll, "vol", dv.pressure) \ - or check_range_local(dcoll, "vol", dv.pressure, .9999999999, 1.00000001): + if check_naninf_local(dcoll, "vol", dv.pressure): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") + if check_range_local(dcoll, "vol", dv.pressure, .9999999999, 1.00000001): + health_error = True + logger.info(f"{rank=}: Pressure range violation.") + from mirgecom.simutil import compare_fluid_solutions component_errors = compare_fluid_solutions(dcoll, state, exact) exittol = .09 @@ -318,7 +330,8 @@ 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, - boundaries=boundaries, gas_model=gas_model) + boundaries=boundaries, gas_model=gas_model, + quadrature_tag=quadrature_tag, use_esdg=use_esdg) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, t_final, constant_cfl) @@ -360,16 +373,26 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") parser.add_argument("--leap", action="store_true", help="use leap timestepper") + parser.add_argument("--esdg", action="store_true", + help="use entropy-stable Euler operator") + parser.add_argument("--overintegration", action="store_true", + help="use over-integration") parser.add_argument("--numpy", action="store_true", help="use numpy-based eager actx.") parser.add_argument("--restart_file", help="root name of restart file") parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -381,7 +404,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_leap=args.leap, - casename=casename, rst_filename=rst_filename) + main(actx_class, use_leap=args.leap, + casename=casename, rst_filename=rst_filename, + use_esdg=args.esdg, use_overintegration=args.esdg or args.overintegration) # vim: foldmethod=marker diff --git a/examples/mixture-mpi.py b/examples/mixture-mpi.py index 65b415a5c..a6a184c73 100644 --- a/examples/mixture-mpi.py +++ b/examples/mixture-mpi.py @@ -71,9 +71,9 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, use_logmgr=True, +def main(actx_class, use_esdg=False, use_leap=False, casename=None, rst_filename=None, - log_dependent=False): + log_dependent=False, use_overintegration=False): """Drive example.""" if casename is None: casename = "mirgecom" @@ -86,7 +86,7 @@ def main(actx_class, use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -142,7 +142,6 @@ def main(actx_class, use_logmgr=True, dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) - use_overintegration = False if use_overintegration: quadrature_tag = DISCR_TAG_QUAD else: @@ -392,7 +391,7 @@ def my_rhs(t, state): return make_obj_array( [euler_operator(dcoll, state=fluid_state, time=t, boundaries=boundaries, gas_model=gas_model, - quadrature_tag=quadrature_tag), + quadrature_tag=quadrature_tag, use_esdg=use_esdg), 0*tseed]) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, @@ -437,8 +436,10 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") + 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("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -450,6 +451,13 @@ def my_rhs(t, state): warn("Automatically turning off DV logging. MIRGE-Com Issue(578)") log_dependent = False + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -458,11 +466,13 @@ def my_rhs(t, state): if args.casename: casename = args.casename rst_filename = None + if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_leap=args.leap, + main(actx_class, use_leap=args.leap, casename=casename, rst_filename=rst_filename, - log_dependent=log_dependent) + use_overintegration=args.overintegration or args.esdg, + use_esdg=args.esdg, log_dependent=log_dependent) # vim: foldmethod=marker diff --git a/examples/multiple-volumes-mpi.py b/examples/multiple-volumes-mpi.py index db5bdd2b9..d595e916b 100644 --- a/examples/multiple-volumes-mpi.py +++ b/examples/multiple-volumes-mpi.py @@ -81,7 +81,7 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, use_logmgr=True, +def main(actx_class, use_esdg=False, use_overintegration=False, use_leap=False, casename=None, rst_filename=None): """Drive the example.""" @@ -96,7 +96,7 @@ def main(actx_class, use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -344,7 +344,7 @@ def my_rhs(t, state): dcoll, state=fluid_state, time=t, boundaries={dd.trace(BTAG_ALL).domain_tag: wall}, gas_model=gas_model, quadrature_tag=quadrature_tag, - dd=dd, comm_tag=dd) + dd=dd, comm_tag=dd, use_esdg=use_esdg) for dd, fluid_state in zip(volume_dds, fluid_states)]) current_dt = my_get_timestep( @@ -384,8 +384,8 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") + parser.add_argument("--esdg", action="store_true", + help="use entropy-stable DG for inviscid terms") parser.add_argument("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -394,6 +394,14 @@ def my_rhs(t, state): parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -405,8 +413,9 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_overintegration=args.overintegration, - use_leap=args.leap, + main(actx_class, + use_leap=args.leap, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker diff --git a/examples/nsmix-mpi.py b/examples/nsmix-mpi.py index fa79a9973..f1bd0fe06 100644 --- a/examples/nsmix-mpi.py +++ b/examples/nsmix-mpi.py @@ -77,14 +77,11 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(use_logmgr=True, +def main(actx_class, use_esdg=False, use_leap=False, casename=None, - rst_filename=None, actx_class=None, + rst_filename=None, log_dependent=True, use_overintegration=False): """Drive example.""" - if actx_class is None: - raise RuntimeError("Array context class missing.") - if casename is None: casename = "mirgecom" @@ -96,7 +93,7 @@ def main(use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -282,7 +279,8 @@ def _ns_operator_for_viz(fluid_state, time): ns_rhs, grad_cv, grad_t = \ ns_operator(dcoll, state=fluid_state, time=time, boundaries=visc_bnds, gas_model=gas_model, - return_gradients=True, quadrature_tag=quadrature_tag) + return_gradients=True, quadrature_tag=quadrature_tag, + use_esdg=use_esdg) return make_obj_array([ns_rhs, grad_cv, grad_t]) get_temperature_update = actx.compile(_get_temperature_update) @@ -552,7 +550,7 @@ def my_rhs(t, state): boundaries=visc_bnds, gas_model=gas_model, gradient_numerical_flux_func=grad_num_flux_func, viscous_numerical_flux_func=viscous_num_flux_func, - quadrature_tag=quadrature_tag) + quadrature_tag=quadrature_tag, use_esdg=use_esdg) cv_rhs = ns_rhs + pyro_eos.get_species_source_terms(cv, fluid_state.temperature) return make_obj_array([cv_rhs, 0*tseed]) @@ -581,7 +579,7 @@ def my_rhs(t, state): ns_rhs, grad_cv, grad_t = \ ns_operator(dcoll, state=current_state, time=current_t, boundaries=visc_bnds, gas_model=gas_model, - return_gradients=True) + return_gradients=True, use_esdg=use_esdg) grad_v = velocity_gradient(current_state.cv, grad_cv) chem_rhs = \ pyro_eos.get_species_source_terms(current_state.cv, @@ -613,10 +611,10 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") parser.add_argument("--leap", action="store_true", help="use leap timestepper") + parser.add_argument("--esdg", action="store_true", + help="use flux-differencing/entropy stable DG for inviscid computations.") parser.add_argument("--numpy", action="store_true", help="use numpy-based eager actx.") parser.add_argument("--restart_file", help="root name of restart file") @@ -624,6 +622,13 @@ def my_rhs(t, state): args = parser.parse_args() from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + warn("Automatically turning off DV logging. MIRGE-Com Issue(578)") log_dependent = False @@ -638,9 +643,9 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(use_logmgr=args.log, use_leap=args.leap, - casename=casename, rst_filename=rst_filename, actx_class=actx_class, - log_dependent=log_dependent, - use_overintegration=args.overintegration) + main(actx_class, use_leap=args.leap, + casename=casename, rst_filename=rst_filename, + log_dependent=log_dependent, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg) # vim: foldmethod=marker diff --git a/examples/poiseuille-local_dt-mpi.py b/examples/poiseuille-local_dt-mpi.py index 3499f7f0e..957aa5d17 100644 --- a/examples/poiseuille-local_dt-mpi.py +++ b/examples/poiseuille-local_dt-mpi.py @@ -71,14 +71,9 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(use_logmgr=True, - use_overintegration=False, - use_leap=False, casename=None, - rst_filename=None, actx_class=None): +def main(actx_class, use_esdg=False, use_overintegration=False, + use_leap=False, casename=None, rst_filename=None): """Drive the example.""" - if actx_class is None: - raise RuntimeError("Array context class missing.") - if casename is None: casename = "mirgecom" @@ -90,7 +85,7 @@ def main(use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -421,7 +416,7 @@ def my_post_step(step, t, dt, state): def my_rhs(t, state): fluid_state = make_fluid_state(state, gas_model) return ns_operator(dcoll, gas_model=gas_model, boundaries=boundaries, - state=fluid_state, time=t, + state=fluid_state, time=t, use_esdg=use_esdg, quadrature_tag=quadrature_tag) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, @@ -478,12 +473,22 @@ def my_rhs(t, state): help="turn on logging") parser.add_argument("--leap", action="store_true", help="use leap timestepper") + parser.add_argument("--esdg", action="store_true", + help="use flux-differencing/entropy stable DG for inviscid computations.") parser.add_argument("--numpy", action="store_true", help="use numpy-based eager actx.") parser.add_argument("--restart_file", help="root name of restart file") parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -495,8 +500,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(use_logmgr=args.log, use_leap=args.leap, - use_overintegration=args.overintegration, - casename=casename, rst_filename=rst_filename, actx_class=actx_class) + main(actx_class, use_leap=args.leap, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, + casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker diff --git a/examples/poiseuille-mpi.py b/examples/poiseuille-mpi.py index 7d850b963..a3838b5b7 100644 --- a/examples/poiseuille-mpi.py +++ b/examples/poiseuille-mpi.py @@ -45,7 +45,8 @@ from mirgecom.steppers import advance_state from mirgecom.boundary import ( PrescribedFluidBoundary, - AdiabaticNoslipWallBoundary + AdiabaticNoslipWallBoundary, + IsothermalWallBoundary ) from mirgecom.transport import SimpleTransport from mirgecom.eos import IdealSingleGas @@ -70,14 +71,9 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(use_logmgr=True, - use_overintegration=False, - use_leap=False, casename=None, - rst_filename=None, actx_class=None): +def main(actx_class, use_esdg=False, use_overintegration=False, + use_leap=False, casename=None, rst_filename=None): """Drive the example.""" - if actx_class is None: - raise RuntimeError("Array context class missing.") - if casename is None: casename = "mirgecom" @@ -89,7 +85,7 @@ def main(use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -109,7 +105,7 @@ def main(use_logmgr=True, # some i/o frequencies nstatus = 1 nviz = 1 - nrestart = 100 + nrestart = 10 nhealth = 1 # some geometry setup @@ -146,8 +142,12 @@ def main(use_logmgr=True, generate_mesh) local_nelements = local_mesh.nelements + # from meshmode.mesh.processing import rotate_mesh_around_axis + # local_mesh = rotate_mesh_around_axis(local_mesh, theta=-np.pi/4) + order = 2 - dcoll = create_discretization_collection(actx, local_mesh, order=order) + dcoll = create_discretization_collection(actx, local_mesh, order=order, + quadrature_order=order+2) nodes = actx.thaw(dcoll.nodes()) if use_overintegration: @@ -198,10 +198,10 @@ def poiseuille_2d(x_vec, eos, cv=None, **kwargs): velocity = make_obj_array([u_x, u_y]) ke = .5*np.dot(velocity, velocity)*mass gamma = eos.gamma() - if cv is not None: - mass = cv.mass - vel = cv.velocity - ke = .5*np.dot(vel, vel)*mass + # if cv is not None: + # mass = cv.mass + # vel = cv.velocity + # ke = .5*np.dot(vel, vel)*mass rho_e = p_x/(gamma-1) + ke return make_conserved(2, mass=mass, energy=rho_e, @@ -210,7 +210,9 @@ def poiseuille_2d(x_vec, eos, cv=None, **kwargs): initializer = poiseuille_2d gas_model = GasModel(eos=IdealSingleGas(), transport=SimpleTransport(viscosity=mu)) - exact = initializer(x_vec=nodes, eos=gas_model.eos) + + from mirgecom.simutil import force_evaluation + exact = force_evaluation(actx, initializer(x_vec=nodes, eos=gas_model.eos)) def _boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context @@ -219,6 +221,12 @@ def _boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, cv=state_minus.cv, **kwargs), gas_model) + use_adiabatic = False + if use_adiabatic: + wall_boundary = AdiabaticNoslipWallBoundary() + else: + wall_boundary = IsothermalWallBoundary(wall_temperature=300) + boundaries = { BoundaryDomainTag("-1"): PrescribedFluidBoundary( @@ -226,8 +234,8 @@ def _boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): BoundaryDomainTag("+1"): PrescribedFluidBoundary( boundary_state_func=_boundary_solution), - BoundaryDomainTag("-2"): AdiabaticNoslipWallBoundary(), - BoundaryDomainTag("+2"): AdiabaticNoslipWallBoundary()} + BoundaryDomainTag("-2"): wall_boundary, + BoundaryDomainTag("+2"): wall_boundary} if rst_filename: current_t = restart_data["t"] @@ -240,7 +248,12 @@ def _boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): # Set the current state from time 0 current_cv = exact - current_state = make_fluid_state(cv=current_cv, gas_model=gas_model) + def _make_fluid_state(cv): + return make_fluid_state(cv=cv, gas_model=gas_model) + + make_fluid_state_compiled = actx.compile(_make_fluid_state) + + current_state = make_fluid_state_compiled(current_cv) vis_timer = None @@ -311,7 +324,7 @@ def my_health_check(state, dv, component_errors): logger.info(f"{rank=}: NANs/Infs in pressure data.") if global_reduce(check_range_local(dcoll, "vol", dv.pressure, 9.999e4, - 1.00101e5), op="lor"): + 1.00101e5), op="lor"): health_error = True from grudge.op import nodal_max, nodal_min p_min = actx.to_numpy(nodal_min(dcoll, "vol", dv.pressure)) @@ -323,14 +336,14 @@ def my_health_check(state, dv, component_errors): logger.info(f"{rank=}: NANs/INFs in temperature data.") if global_reduce(check_range_local(dcoll, "vol", dv.temperature, 348, 350), - op="lor"): + op="lor"): health_error = True from grudge.op import nodal_max, nodal_min t_min = actx.to_numpy(nodal_min(dcoll, "vol", dv.temperature)) t_max = actx.to_numpy(nodal_max(dcoll, "vol", dv.temperature)) logger.info(f"Temperature range violation ({t_min=}, {t_max=})") - exittol = .1 + exittol = 0.1 if max(component_errors) > exittol: health_error = True if rank == 0: @@ -339,8 +352,10 @@ def my_health_check(state, dv, component_errors): return health_error def my_pre_step(step, t, dt, state): - fluid_state = make_fluid_state(cv=state, gas_model=gas_model) + + fluid_state = make_fluid_state_compiled(cv=state) dv = fluid_state.dv + try: component_errors = None @@ -399,9 +414,13 @@ def my_post_step(step, t, dt, state): def my_rhs(t, state): fluid_state = make_fluid_state(state, gas_model) - return ns_operator(dcoll, gas_model=gas_model, boundaries=boundaries, - state=fluid_state, time=t, - quadrature_tag=quadrature_tag) + return ns_operator(dcoll, gas_model=gas_model, + boundaries=boundaries, + state=fluid_state, time=t, use_esdg=use_esdg, + quadrature_tag=quadrature_tag) + + from mirgecom.simutil import force_evaluation + current_state = force_evaluation(actx, current_state) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, t_final, constant_cfl) @@ -417,6 +436,7 @@ def my_rhs(t, state): # Dump the final data if rank == 0: logger.info("Checkpointing final state ...") + final_dv = current_state.dv final_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, t_final, constant_cfl) @@ -445,10 +465,10 @@ def my_rhs(t, state): help="use overintegration in the RHS computations") parser.add_argument("--lazy", action="store_true", help="switch to a lazy computation mode") + parser.add_argument("--esdg", action="store_true", + help="use flux-differencing/entropy stable DG for inviscid computations.") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") parser.add_argument("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -457,6 +477,14 @@ def my_rhs(t, state): parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -468,8 +496,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(use_logmgr=args.log, use_leap=args.leap, - use_overintegration=args.overintegration, - casename=casename, rst_filename=rst_filename, actx_class=actx_class) + main(actx_class, use_leap=args.leap, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, + casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker diff --git a/examples/poiseuille-multispecies-mpi.py b/examples/poiseuille-multispecies-mpi.py new file mode 100644 index 000000000..570b310cb --- /dev/null +++ b/examples/poiseuille-multispecies-mpi.py @@ -0,0 +1,534 @@ +"""Demonstrate a planar Poiseuille flow example with multispecies.""" + +__copyright__ = """ +Copyright (C) 2020 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 logging +import numpy as np +from pytools.obj_array import make_obj_array +from functools import partial + +from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa + +from grudge.eager import EagerDGDiscretization +from grudge.shortcuts import make_visualizer +from grudge.dof_desc import BoundaryDomainTag + +from mirgecom.fluid import make_conserved +from mirgecom.navierstokes import ns_operator +from mirgecom.simutil import get_sim_timestep + +from mirgecom.io import make_init_message +from mirgecom.mpi import mpi_entry_point +from mirgecom.integrators import rk4_step +from mirgecom.steppers import advance_state +from mirgecom.boundary import ( + PrescribedFluidBoundary, + IsothermalWallBoundary +) +from mirgecom.transport import SimpleTransport +from mirgecom.eos import IdealSingleGas # , PyrometheusMixture +from mirgecom.gas_model import GasModel, make_fluid_state +from logpyle import IntervalTimer, set_dt +from mirgecom.euler import extract_vars_for_logging, units_for_logging +from mirgecom.logging_quantities import ( + initialize_logmgr, + logmgr_add_many_discretization_quantities, + logmgr_add_device_name, + logmgr_add_device_memory_usage, + set_sim_state +) + + +logger = logging.getLogger(__name__) + + +class MyRuntimeError(RuntimeError): + """Simple exception to kill the simulation.""" + + pass + + +# Box grid generator widget lifted from @majosm and slightly bent +def _get_box_mesh(dim, a, b, n, t=None): + dim_names = ["x", "y", "z"] + bttf = {} + for i in range(dim): + bttf["-"+str(i+1)] = ["-"+dim_names[i]] + bttf["+"+str(i+1)] = ["+"+dim_names[i]] + from meshmode.mesh.generation import generate_regular_rect_mesh as gen + return gen(a=a, b=b, n=n, boundary_tag_to_face=bttf, mesh_type=t) + + +@mpi_entry_point +def main(actx_class, use_overintegration=False, use_leap=False, casename=None, + rst_filename=None, use_esdg=False): + """Drive the example.""" + if casename is None: + casename = "mirgecom" + + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nparts = comm.Get_size() + + from mirgecom.simutil import global_reduce as _global_reduce + global_reduce = partial(_global_reduce, comm=comm) + + logmgr = initialize_logmgr(True, + filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) + + from mirgecom.array_context import initialize_actx, actx_class_is_profiling + actx = initialize_actx(actx_class, comm) + queue = getattr(actx, "queue", None) + use_profiling = actx_class_is_profiling(actx_class) + + # timestepping control + timestepper = rk4_step + t_final = 1e-7 + current_cfl = 0.05 + current_dt = 1e-8 + current_t = 0 + constant_cfl = False + current_step = 0 + + # some i/o frequencies + nstatus = 100 + nviz = 10 + nrestart = 1000 + nhealth = 100 + + # some geometry setup + dim = 2 + if dim != 2: + raise ValueError("This example must be run with dim = 2.") + x_ch = 1e-4 + left_boundary_location = 0 + right_boundary_location = 0.02 + ybottom = 0. + ytop = .002 + xlen = right_boundary_location - left_boundary_location + ylen = ytop - ybottom + n_refine = 1 + npts_x = n_refine*int(xlen / x_ch) + npts_y = n_refine*int(ylen / x_ch) + + rst_path = "restart_data/" + rst_pattern = ( + rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" + ) + if rst_filename: # read the grid from restart data + rst_filename = f"{rst_filename}-{rank:04d}.pkl" + + from mirgecom.restart import read_restart_data + restart_data = read_restart_data(actx, rst_filename) + local_mesh = restart_data["local_mesh"] + local_nelements = local_mesh.nelements + global_nelements = restart_data["global_nelements"] + assert restart_data["nparts"] == nparts + else: # generate the grid from scratch + npts_axis = (npts_x, npts_y) + box_ll = (left_boundary_location, ybottom) + box_ur = (right_boundary_location, ytop) + generate_mesh = partial(_get_box_mesh, 2, a=box_ll, b=box_ur, n=npts_axis) + from mirgecom.simutil import generate_and_distribute_mesh + 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 = 2 + dcoll = EagerDGDiscretization( + actx, local_mesh, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: default_simplex_group_factory( + base_dim=local_mesh.dim, order=order), + DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order + 1) + }, + mpi_communicator=comm + ) + nodes = actx.thaw(dcoll.nodes()) + + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None + + if logmgr: + logmgr_add_device_name(logmgr, queue) + logmgr_add_device_memory_usage(logmgr, queue) + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, + extract_vars_for_logging, units_for_logging) + + logmgr.add_watches([ + ("step.max", "step = {value}, "), + ("t_sim.max", "sim time: {value:1.6e} s\n"), + ("min_pressure", "------- P (min, max) (Pa) = ({value:1.9e}, "), + ("max_pressure", "{value:1.9e})\n"), + ("min_temperature", "------- T (min, max) (K) = ({value:1.9e}, "), + ("max_temperature", "{value:1.9e})\n"), + ("t_step.max", "------- step walltime: {value:6g} s, "), + ("t_log.max", "log walltime: {value:6g} s") + ]) + + vis_timer = IntervalTimer("t_vis", "Time spent visualizing") + logmgr.add_quantity(vis_timer) + + base_pressure = 100000.0 + pressure_ratio = 1.08 + # MikeA: mu=5e-4, spec_d=1e-4, dt=1e-8, kappa=1e-5 + mu = 5e-4 + kappa = 0. + nspecies = 2 + species_diffusivity = 1e-5 * np.ones(nspecies) + xlen = right_boundary_location - left_boundary_location + ylen = ytop - ybottom + + def poiseuille_2d(x_vec, eos, cv=None, **kwargs): + y = x_vec[1] + x = x_vec[0] + # zeros = 0*x + ones = 0*x + 1. + x0 = left_boundary_location + xmax = right_boundary_location + xlen = xmax - x0 + wgt1 = actx.np.less(x, xlen/2) + wgt2 = 1 - wgt1 + # xcor = x*ones + # leno2 = xlen/2*ones + p_low = base_pressure + p_hi = pressure_ratio*base_pressure + dp = p_hi - p_low + dpdx = dp/xlen + h = ytop - ybottom + u_x = dpdx*y*(h - y)/(2*mu) + print(f"flow speed = {dpdx*h*h/(8*mu)}") + p_x = p_hi - dpdx*x + rho = 1.0*ones + mass = 0*x + rho + u_y = 0*x + velocity = make_obj_array([u_x, u_y]) + ke = .5*np.dot(velocity, velocity)*mass + gamma = eos.gamma() + rho_y = rho * make_obj_array([1.0/nspecies for _ in range(nspecies)]) + rho_y[0] = wgt1*rho_y[0] + rho_y[1] = wgt2*rho_y[1] + if cv is not None: + rho_y = wgt1*rho_y + wgt2*mass*cv.species_mass_fractions + + rho_e = p_x/(gamma-1) + ke + return make_conserved(2, mass=mass, energy=rho_e, + momentum=mass*velocity, + species_mass=rho_y) + + initializer = poiseuille_2d + gas_model = GasModel(eos=IdealSingleGas(), + transport=SimpleTransport( + viscosity=mu, thermal_conductivity=kappa, + species_diffusivity=species_diffusivity)) + exact = initializer(x_vec=nodes, eos=gas_model.eos) + + def _exact_boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): + actx = state_minus.array_context + bnd_discr = dcoll.discr_from_dd(dd_bdry) + nodes = actx.thaw(bnd_discr.nodes()) + return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, + cv=state_minus.cv, **kwargs), gas_model) + + boundaries = { + BoundaryDomainTag("-1"): + PrescribedFluidBoundary(boundary_state_func=_exact_boundary_solution), + BoundaryDomainTag("+1"): + PrescribedFluidBoundary(boundary_state_func=_exact_boundary_solution), + BoundaryDomainTag("-2"): + IsothermalWallBoundary(wall_temperature=348.5), + BoundaryDomainTag("+2"): + IsothermalWallBoundary(wall_temperature=348.5)} + + if rst_filename: + current_t = restart_data["t"] + current_step = restart_data["step"] + current_cv = restart_data["cv"] + if logmgr: + from mirgecom.logging_quantities import logmgr_set_time + logmgr_set_time(logmgr, current_step, current_t) + else: + # Set the current state from time 0 + current_cv = exact + + vis_timer = None + + visualizer = make_visualizer(dcoll, order) + + eosname = gas_model.eos.__class__.__name__ + init_message = make_init_message(dim=dim, order=order, + nelements=local_nelements, + global_nelements=global_nelements, + dt=current_dt, t_final=t_final, nstatus=nstatus, + nviz=nviz, cfl=current_cfl, + constant_cfl=constant_cfl, initname=casename, + eosname=eosname, casename=casename) + if rank == 0: + logger.info(init_message) + + def my_write_status(step, t, dt, state, component_errors): + dv = state.dv + from grudge.op import nodal_min, nodal_max + p_min = actx.to_numpy(nodal_min(dcoll, "vol", dv.pressure)) + p_max = actx.to_numpy(nodal_max(dcoll, "vol", dv.pressure)) + t_min = actx.to_numpy(nodal_min(dcoll, "vol", dv.temperature)) + t_max = actx.to_numpy(nodal_max(dcoll, "vol", dv.temperature)) + if constant_cfl: + cfl = current_cfl + else: + from mirgecom.viscous import get_viscous_cfl + cfl = actx.to_numpy(nodal_max(dcoll, "vol", + get_viscous_cfl(dcoll, dt, state))) + if rank == 0: + logger.info(f"Step: {step}, T: {t}, DT: {dt}, CFL: {cfl}\n" + f"----- Pressure({p_min}, {p_max})\n" + f"----- Temperature({t_min}, {t_max})\n" + "----- errors=" + + ", ".join("%.3g" % en for en in component_errors)) + + def my_write_viz(step, t, state, dv): + resid = state - exact + viz_fields = [("cv", state), + ("dv", dv), + ("poiseuille", exact), + ("resid", resid)] + + from mirgecom.simutil import write_visfile + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, + step=step, t=t, overwrite=True) + + def my_write_restart(step, t, state): + rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "cv": state, + "t": t, + "step": step, + "order": order, + "global_nelements": global_nelements, + "num_parts": nparts + } + from mirgecom.restart import write_restart_file + write_restart_file(actx, rst_data, rst_fname, comm) + + def my_health_check(state, dv, component_errors): + health_error = False + from mirgecom.simutil import check_naninf_local, check_range_local + if check_naninf_local(dcoll, "vol", dv.pressure): + health_error = True + logger.info(f"{rank=}: NANs/Infs in pressure data.") + + if global_reduce(check_range_local(dcoll, "vol", dv.pressure, 9.999e4, + 1.00101e5), op="lor"): + health_error = False + from grudge.op import nodal_max, nodal_min + p_min = actx.to_numpy(nodal_min(dcoll, "vol", dv.pressure)) + p_max = actx.to_numpy(nodal_max(dcoll, "vol", dv.pressure)) + logger.info(f"Pressure range violation ({p_min=}, {p_max=})") + + if check_naninf_local(dcoll, "vol", dv.temperature): + health_error = True + logger.info(f"{rank=}: NANs/INFs in temperature data.") + + if global_reduce(check_range_local(dcoll, "vol", dv.temperature, 348, 350), + op="lor"): + health_error = False + from grudge.op import nodal_max, nodal_min + t_min = actx.to_numpy(nodal_min(dcoll, "vol", dv.temperature)) + t_max = actx.to_numpy(nodal_max(dcoll, "vol", dv.temperature)) + logger.info(f"Temperature range violation ({t_min=}, {t_max=})") + + exittol = 1e7 + if max(component_errors) > exittol: + health_error = False + if rank == 0: + logger.info("Solution diverged from exact soln.") + + return health_error + + def my_pre_step(step, t, dt, state): + fluid_state = make_fluid_state(cv=state, gas_model=gas_model) + dv = fluid_state.dv + try: + component_errors = None + + if logmgr: + logmgr.tick_before() + + from mirgecom.simutil import check_step + do_viz = check_step(step=step, interval=nviz) + do_restart = check_step(step=step, interval=nrestart) + do_health = check_step(step=step, interval=nhealth) + do_status = check_step(step=step, interval=nstatus) + + if do_health: + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(dcoll, state, exact) + health_errors = global_reduce( + my_health_check(state, dv, component_errors), op="lor") + if health_errors: + if rank == 0: + logger.info("Fluid solution failed health check.") + raise MyRuntimeError("Failed simulation health check.") + + if do_restart: + my_write_restart(step=step, t=t, state=state) + + if do_viz: + my_write_viz(step=step, t=t, state=state, dv=dv) + + dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, + t_final, constant_cfl) + + if do_status: # needed because logging fails to make output + if component_errors is None: + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(dcoll, state, exact) + my_write_status(step=step, t=t, dt=dt, state=fluid_state, + component_errors=component_errors) + + except MyRuntimeError: + if rank == 0: + logger.info("Errors detected; attempting graceful exit.") + my_write_viz(step=step, t=t, state=state, dv=dv) + my_write_restart(step=step, t=t, state=state) + raise + + return state, dt + + def my_post_step(step, t, dt, state): + # Logmgr needs to know about EOS, dt, dim? + # imo this is a design/scope flaw + if logmgr: + set_dt(logmgr, dt) + set_sim_state(logmgr, dim, state, gas_model.eos) + logmgr.tick_after() + return state, dt + + orig = np.zeros(shape=(dim,)) + orig[0] = 2*xlen/5. + orig[1] = 7*ylen/10. + + def acoustic_pulse(time, fluid_cv, gas_model): + from mirgecom.initializers import AcousticPulse + acoustic_pulse = AcousticPulse(dim=dim, amplitude=5000.0, width=.0001, + center=orig) + # return fluid_cv + return acoustic_pulse(nodes, cv=fluid_cv, eos=gas_model.eos) + + def my_rhs(t, state): + fluid_state = make_fluid_state(state, gas_model) + return ns_operator(dcoll, gas_model=gas_model, boundaries=boundaries, + state=fluid_state, time=t, use_esdg=use_esdg, + quadrature_tag=quadrature_tag) + + current_state = make_fluid_state( + cv=acoustic_pulse(current_t, current_cv, gas_model), gas_model=gas_model) + + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, + current_cfl, t_final, constant_cfl) + + current_step, current_t, current_cv = \ + advance_state(rhs=my_rhs, timestepper=timestepper, + pre_step_callback=my_pre_step, + post_step_callback=my_post_step, dt=current_dt, + state=current_state.cv, t=current_t, t_final=t_final) + + current_state = make_fluid_state(cv=current_cv, gas_model=gas_model) + + # Dump the final data + if rank == 0: + logger.info("Checkpointing final state ...") + final_dv = current_state.dv + final_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, + current_cfl, t_final, constant_cfl) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(dcoll, current_state.cv, exact) + + my_write_viz(step=current_step, t=current_t, state=current_state.cv, dv=final_dv) + my_write_restart(step=current_step, t=current_t, state=current_state) + my_write_status(step=current_step, t=current_t, dt=final_dt, + state=current_state, component_errors=component_errors) + + if logmgr: + logmgr.close() + elif use_profiling: + print(actx.tabulate_profiling_data()) + + finish_tol = 1e-16 + assert np.abs(current_t - t_final) < finish_tol + + +if __name__ == "__main__": + import argparse + casename = "poiseuille" + 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("--lazy", action="store_true", + help="switch to a lazy computation mode") + parser.add_argument("--profiling", action="store_true", + help="turn on detailed performance profiling") + parser.add_argument("--leap", action="store_true", + help="use leap timestepper") + parser.add_argument("--esdg", action="store_true", + help="use flux-differencing/entropy stable DG for inviscid computations.") + parser.add_argument("--numpy", action="store_true", + help="use numpy-based eager actx.") + parser.add_argument("--restart_file", help="root name of restart file") + parser.add_argument("--casename", help="casename to use for i/o") + args = parser.parse_args() + + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + + from mirgecom.array_context import get_reasonable_array_context_class + actx_class = get_reasonable_array_context_class( + lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) + + logging.basicConfig(format="%(message)s", level=logging.INFO) + if args.casename: + casename = args.casename + rst_filename = None + if args.restart_file: + rst_filename = args.restart_file + + main(actx_class, use_leap=args.leap, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, + casename=casename, rst_filename=rst_filename) + +# vim: foldmethod=marker diff --git a/examples/pulse-mpi.py b/examples/pulse-mpi.py index 7953e005a..0337cbc2c 100644 --- a/examples/pulse-mpi.py +++ b/examples/pulse-mpi.py @@ -77,7 +77,7 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, use_logmgr=True, +def main(actx_class, use_esdg=False, use_overintegration=False, use_leap=False, casename=None, rst_filename=None): """Drive the example.""" @@ -92,7 +92,7 @@ def main(actx_class, use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -329,7 +329,7 @@ def my_rhs(t, state): fluid_state = make_fluid_state(cv=state, gas_model=gas_model) return euler_operator(dcoll, state=fluid_state, time=t, boundaries=boundaries, - gas_model=gas_model, + gas_model=gas_model, use_esdg=use_esdg, quadrature_tag=quadrature_tag) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, @@ -369,16 +369,24 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") parser.add_argument("--leap", action="store_true", help="use leap timestepper") + parser.add_argument("--esdg", action="store_true", + help="use entropy-stable dg for inviscid terms.") parser.add_argument("--numpy", action="store_true", help="use numpy-based eager actx.") parser.add_argument("--restart_file", help="root name of restart file") parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -390,7 +398,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_overintegration=args.overintegration, + main(actx_class, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, use_leap=args.leap, casename=casename, rst_filename=rst_filename) diff --git a/examples/scalar-advdiff-mpi.py b/examples/scalar-advdiff-mpi.py index b18c260d8..75ddb1f77 100644 --- a/examples/scalar-advdiff-mpi.py +++ b/examples/scalar-advdiff-mpi.py @@ -36,10 +36,12 @@ from mirgecom.transport import SimpleTransport from mirgecom.navierstokes import ns_operator from mirgecom.simutil import ( - get_sim_timestep, + # get_sim_timestep, generate_and_distribute_mesh, compare_fluid_solutions ) +from mirgecom.limiter import bound_preserving_limiter +from mirgecom.fluid import make_conserved from mirgecom.io import make_init_message from mirgecom.mpi import mpi_entry_point @@ -69,9 +71,8 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, use_logmgr=True, - use_leap=False, casename=None, - rst_filename=None): +def main(actx_class, use_overintegration=False, use_esdg=False, + use_leap=False, casename=None, rst_filename=None): """Drive example.""" if casename is None: casename = "mirgecom" @@ -84,7 +85,7 @@ def main(actx_class, use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -99,21 +100,22 @@ def main(actx_class, use_logmgr=True, timestepper = RK4MethodBuilder("state") else: timestepper = rk4_step - t_final = 0.1 + + t_final = 2e-4 current_cfl = 0.1 - current_dt = .001 + current_dt = 1e-5 current_t = 0 - constant_cfl = True + constant_cfl = False # some i/o frequencies - nstatus = 1 - nrestart = 5 - nviz = 100 - nhealth = 1 + nstatus = 100 + nrestart = 100 + nviz = 10 + nhealth = 100 dim = 2 - nel_1d = 4 - order = 1 + nel_1d = 8 + order = 3 rst_path = "restart_data/" rst_pattern = ( @@ -142,6 +144,47 @@ def main(actx_class, use_logmgr=True, dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) + from grudge.dof_desc import DISCR_TAG_QUAD + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None + + def _limit_fluid_cv(cv, pressure, temperature, dd=None): + # if True: + # return cv + actx = cv.array_context + # limit species + spec_lim = make_obj_array([ + bound_preserving_limiter(dcoll, cv.species_mass_fractions[i], mmin=0.0, + dd=dd) + for i in range(nspecies) + ]) + spec_lim = actx.np.where(actx.np.greater(spec_lim, 0.0), spec_lim, 0.0) + + # normalize to ensure sum_Yi = 1.0 + # aux = cv.mass*0.0 + # for i in range(0, nspecies): + # aux = aux + spec_lim[i] + # spec_lim = spec_lim/aux + + # recompute density + # mass_lim = eos.get_density(pressure=pressure, + # temperature=temperature, species_mass_fractions=spec_lim) + + # recompute energy + # energy_lim = mass_lim*(gas_model.eos.get_internal_energy( + # temperature, species_mass_fractions=spec_lim) + # + 0.5*np.dot(cv.velocity, cv.velocity) + # ) + + # make a new CV with the limited variables + return make_conserved(dim=dim, mass=cv.mass, energy=cv.energy, + momentum=cv.momentum, + species_mass=cv.mass*spec_lim) + use_limiter = False + limiter_function = _limit_fluid_cv if use_limiter else None + def vol_min(x): from grudge.op import nodal_min return actx.to_numpy(nodal_min(dcoll, "vol", x))[()] @@ -176,23 +219,24 @@ def vol_max(x): ]) # soln setup and init - nspecies = 1 + nspecies = 4 centers = make_obj_array([np.zeros(shape=(dim,)) for i in range(nspecies)]) velocity = np.zeros(shape=(dim,)) - velocity[0] = 1. + velocity[0] = 300. wave_vector = np.zeros(shape=(dim,)) wave_vector[0] = 1. wave_vector = wave_vector / np.sqrt(np.dot(wave_vector, wave_vector)) - spec_y0s = np.ones(shape=(nspecies,)) + spec_y0s = np.zeros(shape=(nspecies,)) spec_amplitudes = np.ones(shape=(nspecies,)) spec_omegas = 2. * np.pi * np.ones(shape=(nspecies,)) - kappa = 1e-5 - sigma = 1e-5 - spec_diff = .1 - spec_diffusivities = spec_diff * np.ones(nspecies) - transport_model = SimpleTransport(viscosity=sigma, thermal_conductivity=kappa, + kappa = 0.0 + mu = 1e-5 + spec_diff = 1e-1 + spec_diffusivities = np.array([spec_diff * 1./float(j+1) + for j in range(nspecies)]) + transport_model = SimpleTransport(viscosity=mu, thermal_conductivity=kappa, species_diffusivity=spec_diffusivities) eos = IdealSingleGas() @@ -201,19 +245,23 @@ def vol_max(x): from mirgecom.initializers import MulticomponentTrig initializer = MulticomponentTrig(dim=dim, nspecies=nspecies, + p0=101325, rho0=1.3, spec_centers=centers, velocity=velocity, spec_y0s=spec_y0s, spec_amplitudes=spec_amplitudes, spec_omegas=spec_omegas, spec_diffusivities=spec_diffusivities, - wave_vector=wave_vector) + wave_vector=wave_vector, + trig_function=actx.np.sin) def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): actx = state_minus.array_context bnd_discr = dcoll.discr_from_dd(dd_bdry) nodes = actx.thaw(bnd_discr.nodes()) return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, - **kwargs), gas_model) + **kwargs), gas_model, + limiter_func=limiter_function, + limiter_dd=dd_bdry) boundaries = {} @@ -228,7 +276,8 @@ def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): # Set the current state from time 0 current_cv = initializer(nodes) - current_state = make_fluid_state(current_cv, gas_model) + current_state = make_fluid_state(current_cv, gas_model, + limiter_func=limiter_function) convective_speed = np.sqrt(np.dot(velocity, velocity)) c = current_state.speed_of_sound mach = vol_max(convective_speed / c) @@ -285,11 +334,14 @@ def my_write_restart(step, t, cv): 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, .99999999, 1.00000001): + if check_naninf_local(dcoll, "vol", pressure): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") + if check_range_local(dcoll, "vol", pressure, 101324.99, 101325.01): + health_error = True + logger.info(f"{rank=}: Pressure out of expected range.") + exittol = .09 if max(component_errors) > exittol: health_error = False @@ -342,8 +394,9 @@ def my_pre_step(step, t, dt, state): logger.info("Errors detected; attempting graceful exit.") raise - dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, t_final, - constant_cfl) + # dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, t_final, + # constant_cfl) + return state, dt def my_post_step(step, t, dt, state): @@ -356,12 +409,15 @@ def my_post_step(step, t, dt, state): return state, dt def my_rhs(t, state): - fluid_state = make_fluid_state(state, gas_model) + fluid_state = make_fluid_state(state, gas_model, + limiter_func=limiter_function) return ns_operator(dcoll, state=fluid_state, time=t, - boundaries=boundaries, gas_model=gas_model) + boundaries=boundaries, gas_model=gas_model, + quadrature_tag=quadrature_tag, use_esdg=use_esdg, + limiter_func=limiter_function) - current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, - current_cfl, t_final, constant_cfl) + # current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, + # current_cfl, t_final, constant_cfl) current_step, current_t, current_cv = \ advance_state(rhs=my_rhs, timestepper=timestepper, @@ -399,16 +455,26 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") parser.add_argument("--leap", action="store_true", help="use leap timestepper") + parser.add_argument("--esdg", action="store_true", + help="use entropy-stable rhs operator") + parser.add_argument("--overintegration", action="store_true", + help="use overintegration") parser.add_argument("--numpy", action="store_true", help="use numpy-based eager actx.") parser.add_argument("--restart_file", help="root name of restart file") parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class(lazy=args.lazy, distributed=True, profiling=args.profiling, @@ -421,7 +487,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_leap=args.leap, + main(actx_class, use_leap=args.leap, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker diff --git a/examples/scalar-lump-mpi.py b/examples/scalar-lump-mpi.py index 912c0cf3a..4a02af72f 100644 --- a/examples/scalar-lump-mpi.py +++ b/examples/scalar-lump-mpi.py @@ -28,6 +28,7 @@ from functools import partial from pytools.obj_array import make_obj_array +from grudge.dof_desc import DISCR_TAG_QUAD from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.shortcuts import make_visualizer @@ -66,8 +67,8 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, use_logmgr=True, - use_leap=False, casename=None, +def main(actx_class, use_leap=False, casename=None, + use_overintegration=False, use_esdg=False, rst_filename=None): """Drive example.""" if casename is None: @@ -81,7 +82,7 @@ def main(actx_class, use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -96,17 +97,18 @@ def main(actx_class, use_logmgr=True, timestepper = RK4MethodBuilder("state") else: timestepper = rk4_step - t_final = 0.005 + + t_final = 2e-2 current_cfl = 1.0 - current_dt = .001 + current_dt = 1e-3 current_t = 0 constant_cfl = False # some i/o frequencies - nstatus = 1 - nrestart = 5 - nviz = 1 - nhealth = 1 + nstatus = 100 + nrestart = 10000 + nviz = 10 + nhealth = 100 dim = 2 rst_path = "restart_data/" @@ -136,6 +138,11 @@ def main(actx_class, use_logmgr=True, dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None + vis_timer = None if logmgr: @@ -247,14 +254,17 @@ def my_write_restart(step, t, state): 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, .99999999, 1.00000001): + if check_naninf_local(dcoll, "vol", pressure): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") + if check_range_local(dcoll, "vol", pressure, .99999999, 1.00000001): + health_error = True + logger.info(f"{rank=}: Pressure range violation.") + exittol = .09 if max(component_errors) > exittol: - health_error = True + # health_error = True if rank == 0: logger.info("Solution diverged from exact soln.") @@ -330,7 +340,8 @@ 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, - boundaries=boundaries, gas_model=gas_model) + boundaries=boundaries, gas_model=gas_model, + quadrature_tag=quadrature_tag, use_esdg=use_esdg) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, t_final, constant_cfl) @@ -372,16 +383,26 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") parser.add_argument("--leap", action="store_true", help="use leap timestepper") + parser.add_argument("--esdg", action="store_true", + help="use entropy-stable Euler operator") + parser.add_argument("--overintegration", action="store_true", + help="use over-integration") parser.add_argument("--numpy", action="store_true", help="use numpy-based eager actx.") parser.add_argument("--restart_file", help="root name of restart file") parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -393,7 +414,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_leap=args.leap, - casename=casename, rst_filename=rst_filename) + main(actx_class, use_leap=args.leap, + casename=casename, rst_filename=rst_filename, + use_esdg=args.esdg, use_overintegration=args.esdg or args.overintegration) # vim: foldmethod=marker diff --git a/examples/sod-mpi.py b/examples/sod-mpi.py index 34902eb9e..e23d22b57 100644 --- a/examples/sod-mpi.py +++ b/examples/sod-mpi.py @@ -65,9 +65,8 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, use_logmgr=True, - use_leap=False, casename=None, - rst_filename=None): +def main(actx_class, use_overintegration=False, use_esdg=False, + use_leap=False, casename=None, rst_filename=None): """Drive the example.""" if casename is None: casename = "mirgecom" @@ -80,7 +79,7 @@ def main(actx_class, use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -94,18 +93,18 @@ def main(actx_class, use_logmgr=True, timestepper = RK4MethodBuilder("state") else: timestepper = rk4_step - t_final = 0.01 - current_cfl = 1.0 - current_dt = .0001 + t_final = 1e-4 + current_cfl = 0.01 + current_dt = 1e-6 current_t = 0 constant_cfl = False current_step = 0 # some i/o frequencies - nstatus = 10 - nrestart = 5 - nviz = 10 - nhealth = 10 + nstatus = 1 + nrestart = 1000 + nviz = 1 + nhealth = 1 dim = 1 rst_path = "restart_data/" @@ -122,19 +121,31 @@ def main(actx_class, use_logmgr=True, assert restart_data["num_parts"] == num_parts else: # generate the grid from scratch from meshmode.mesh.generation import generate_regular_rect_mesh - nel_1d = 24 - box_ll = -5.0 - box_ur = 5.0 + nel_1d = 32 + box_ll = 0.0 + box_ur = 1.0 generate_mesh = partial(generate_regular_rect_mesh, a=(box_ll,)*dim, b=(box_ur,) * dim, nelements_per_axis=(nel_1d,)*dim) local_mesh, global_nelements = generate_and_distribute_mesh(comm, generate_mesh) local_nelements = local_mesh.nelements - order = 1 - dcoll = create_discretization_collection(actx, local_mesh, order=order) + order = 4 + dcoll = create_discretization_collection(actx, local_mesh, order=order, + quadrature_order=order+2) nodes = actx.thaw(dcoll.nodes()) + # TODO: Fix this wonky dt estimate + from grudge.dt_utils import h_min_from_volume + cn = 0.5*(order + 1)**2 + current_dt = current_cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn + + from grudge.dof_desc import DISCR_TAG_QUAD + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None + vis_timer = None if logmgr: @@ -156,7 +167,7 @@ def main(actx_class, use_logmgr=True, ]) initializer = SodShock1D(dim=dim) - eos = IdealSingleGas() + eos = IdealSingleGas(gas_const=1.0) gas_model = GasModel(eos=eos) def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): @@ -169,6 +180,7 @@ def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): boundaries = { BTAG_ALL: PrescribedFluidBoundary(boundary_state_func=boundary_solution) } + if rst_filename: current_t = restart_data["t"] current_step = restart_data["step"] @@ -193,6 +205,7 @@ def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): nviz=nviz, cfl=current_cfl, constant_cfl=constant_cfl, initname=initname, eosname=eosname, casename=casename) + if rank == 0: logger.info(init_message) @@ -211,9 +224,11 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): if resid is None: resid = state - exact viz_fields = [("cv", state), + ("vel", state.velocity), ("dv", dv), ("exact", exact), ("residual", resid)] + from mirgecom.simutil import write_visfile write_visfile(dcoll, viz_fields, visualizer, vizname=casename, step=step, t=t, overwrite=True, vis_timer=vis_timer, @@ -236,17 +251,18 @@ def my_write_restart(state, step, t): 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, .09, 1.1): + # from mirgecom.simutil import check_naninf_local, check_range_local + from mirgecom.simutil import check_naninf_local + # or check_range_local(dcoll, "vol", pressure, .09, 1.1): + if check_naninf_local(dcoll, "vol", pressure): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") - exittol = .09 - if max(component_errors) > exittol: - health_error = True - if rank == 0: - logger.info("Solution diverged from exact soln.") + # exittol = .09 + # if max(component_errors) > exittol: + # health_error = True + # if rank == 0: + # logger.info("Solution diverged from exact soln.") return health_error @@ -320,7 +336,9 @@ def my_post_step(step, t, dt, state): def my_rhs(t, state): fluid_state = make_fluid_state(cv=state, gas_model=gas_model) return euler_operator(dcoll, state=fluid_state, time=t, - boundaries=boundaries, gas_model=gas_model) + boundaries=boundaries, gas_model=gas_model, + quadrature_tag=quadrature_tag, + use_esdg=use_esdg) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, t_final, constant_cfl) @@ -360,16 +378,26 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") parser.add_argument("--leap", action="store_true", help="use leap timestepper") + 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("--numpy", action="store_true", help="use numpy-based eager actx.") parser.add_argument("--restart_file", help="root name of restart file") parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -381,7 +409,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_leap=args.leap, - casename=casename, rst_filename=rst_filename) + main(actx_class, use_leap=args.leap, + casename=casename, rst_filename=rst_filename, + use_overintegration=args.overintegration or args.esdg, use_esdg=args.esdg) # vim: foldmethod=marker diff --git a/examples/taylor-green-mpi.py b/examples/taylor-green-mpi.py new file mode 100644 index 000000000..ff7e290f1 --- /dev/null +++ b/examples/taylor-green-mpi.py @@ -0,0 +1,363 @@ +"""Demonstrate the inviscid Taylor-Green vortex problem.""" + +__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 logging + +from mirgecom.mpi import mpi_entry_point + +from functools import partial + +from grudge.shortcuts import make_visualizer + +from mirgecom.euler import euler_operator +from mirgecom.simutil import ( + generate_and_distribute_mesh +) +from mirgecom.io import make_init_message +from mirgecom.discretization import create_discretization_collection +from mirgecom.integrators import lsrk54_step +from mirgecom.steppers import advance_state +from mirgecom.initializers import InviscidTaylorGreenVortex +from mirgecom.eos import IdealSingleGas +from mirgecom.gas_model import ( + GasModel, + make_fluid_state +) +from logpyle import IntervalTimer, set_dt +from mirgecom.euler import extract_vars_for_logging, units_for_logging +from mirgecom.logging_quantities import ( + initialize_logmgr, + logmgr_add_many_discretization_quantities, + logmgr_add_device_name, + logmgr_add_device_memory_usage, + set_sim_state +) + +logger = logging.getLogger(__name__) + + +class MyRuntimeError(RuntimeError): + """Simple exception to kill the simulation.""" + + pass + + +@mpi_entry_point +def main(actx_class, order=1, t_final=1, resolution=4, + use_overintegration=False, use_esdg=False, + casename=None, rst_filename=None): + """Drive the example.""" + if casename is None: + casename = "mirgecom" + + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + num_parts = comm.Get_size() + + from mirgecom.simutil import global_reduce as _global_reduce + global_reduce = partial(_global_reduce, comm=comm) + + logmgr = initialize_logmgr(True, + filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) + + from mirgecom.array_context import initialize_actx, actx_class_is_profiling + actx = initialize_actx(actx_class, comm) + queue = getattr(actx, "queue", None) + use_profiling = actx_class_is_profiling(actx_class) + + # timestepping control + current_step = 0 + timestepper = lsrk54_step + t_final = 5e-3 + current_cfl = 1.0 + current_dt = 1e-3 + current_t = 0 + constant_cfl = False + + # some i/o frequencies + nstatus = 100000 + nrestart = 100 + nviz = 100 + nhealth = 100 + + dim = 3 + rst_path = "restart_data/" + rst_pattern = ( + rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" + ) + if rst_filename: # read the grid from restart data + rst_filename = f"{rst_filename}-{rank:04d}.pkl" + from mirgecom.restart import read_restart_data + restart_data = read_restart_data(actx, rst_filename) + local_mesh = restart_data["local_mesh"] + local_nelements = local_mesh.nelements + global_nelements = restart_data["global_nelements"] + assert restart_data["num_parts"] == num_parts + else: # generate the grid from scratch + from meshmode.mesh.generation import generate_regular_rect_mesh + box_ll = -np.pi + box_ur = np.pi + generate_mesh = partial(generate_regular_rect_mesh, + a=(box_ll,)*dim, + b=(box_ur,) * dim, + nelements_per_axis=(resolution,)*dim, + periodic=(True,)*dim) + local_mesh, global_nelements = generate_and_distribute_mesh(comm, + generate_mesh) + local_nelements = local_mesh.nelements + + from grudge.dof_desc import DISCR_TAG_QUAD + + dcoll = create_discretization_collection(actx, local_mesh, order=order) + nodes = actx.thaw(dcoll.nodes()) + + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None + + vis_timer = None + + if logmgr: + logmgr_add_device_name(logmgr, queue) + logmgr_add_device_memory_usage(logmgr, queue) + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, + extract_vars_for_logging, units_for_logging) + + vis_timer = IntervalTimer("t_vis", "Time spent visualizing") + logmgr.add_quantity(vis_timer) + + logmgr.add_watches([ + ("step.max", "step = {value}, "), + ("t_sim.max", "sim time: {value:1.6e} s\n"), + ("min_pressure", "------- P (min, max) (Pa) = ({value:1.9e}, "), + ("max_pressure", "{value:1.9e})\n"), + ("t_step.max", "------- step walltime: {value:6g} s, "), + ("t_log.max", "log walltime: {value:6g} s") + ]) + + eos = IdealSingleGas() + gas_model = GasModel(eos=eos) + + # Periodic domain, no boundary conditions (yay!) + boundaries = {} + + initial_condition = InviscidTaylorGreenVortex() + + if rst_filename: + current_t = restart_data["t"] + current_step = restart_data["step"] + current_cv = restart_data["cv"] + if logmgr: + from mirgecom.logging_quantities import logmgr_set_time + logmgr_set_time(logmgr, current_step, current_t) + else: + # Set the current state from time 0 + from mirgecom.simutil import force_evaluation + current_cv = force_evaluation(actx, + initial_condition(x_vec=nodes, eos=eos)) + + visualizer = make_visualizer(dcoll) + + initname = "taylorgreen" + eosname = eos.__class__.__name__ + init_message = make_init_message(dim=dim, order=order, + nelements=local_nelements, + global_nelements=global_nelements, + dt=current_dt, t_final=t_final, nstatus=nstatus, + nviz=nviz, cfl=current_cfl, + constant_cfl=constant_cfl, initname=initname, + eosname=eosname, casename=casename) + if rank == 0: + logger.info(init_message) + + def _make_fluid_state(cv, tseed): + return make_fluid_state(cv=cv, gas_model=gas_model) + + make_fluid_state_compiled = actx.compile(_make_fluid_state) + # current_state = make_fluid_state_compiled(current_cv, gas_model) + + def my_write_viz(step, t, state, dv=None): + if dv is None: + dv = eos.dependent_vars(state) + viz_fields = [("cv", state), + ("dv", dv)] + from mirgecom.simutil import write_visfile + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, + step=step, t=t, overwrite=True, vis_timer=vis_timer) + + def my_write_restart(step, t, state): + rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "cv": state, + "t": t, + "step": step, + "order": order, + "global_nelements": global_nelements, + "num_parts": num_parts + } + from mirgecom.restart import write_restart_file + write_restart_file(actx, rst_data, rst_fname, comm) + + def my_health_check(pressure): + health_error = False + from mirgecom.simutil import check_naninf_local + if check_naninf_local(dcoll, "vol", pressure): + health_error = True + logger.info(f"{rank=}: Invalid pressure data found.") + return health_error + + def my_pre_step(step, t, dt, state): + fluid_state = make_fluid_state_compiled(state, gas_model) + dv = fluid_state.dv + + try: + + if logmgr: + logmgr.tick_before() + + from mirgecom.simutil import check_step + do_viz = check_step(step=step, interval=nviz) + do_restart = check_step(step=step, interval=nrestart) + do_health = check_step(step=step, interval=nhealth) + + if do_health: + health_errors = global_reduce(my_health_check(dv.pressure), op="lor") + if health_errors: + if rank == 0: + logger.info("Fluid solution failed health check.") + raise MyRuntimeError("Failed simulation health check.") + + if do_restart: + my_write_restart(step=step, t=t, state=state) + + if do_viz: + my_write_viz(step=step, t=t, state=state, dv=dv) + + except MyRuntimeError: + if rank == 0: + logger.info("Errors detected; attempting graceful exit.") + my_write_viz(step=step, t=t, state=state) + my_write_restart(step=step, t=t, state=state) + raise + + # dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, t_final, + # constant_cfl) + return state, dt + + def my_post_step(step, t, dt, state): + if logmgr: + set_dt(logmgr, dt) + set_sim_state(logmgr, dim, state, eos) + logmgr.tick_after() + return state, dt + + def my_rhs(t, state): + fluid_state = make_fluid_state(cv=state, gas_model=gas_model) + return euler_operator(dcoll, state=fluid_state, time=t, + boundaries=boundaries, + gas_model=gas_model, use_esdg=use_esdg, + quadrature_tag=quadrature_tag) + + current_step, current_t, current_cv = \ + advance_state(rhs=my_rhs, timestepper=timestepper, + pre_step_callback=my_pre_step, + post_step_callback=my_post_step, dt=current_dt, + state=current_cv, t=current_t, t_final=t_final) + + # Dump the final data + if rank == 0: + logger.info("Checkpointing final state ...") + final_state = make_fluid_state(current_cv, gas_model) + final_dv = final_state.dv + + my_write_viz(step=current_step, t=current_t, state=current_cv, dv=final_dv) + my_write_restart(step=current_step, t=current_t, state=current_cv) + + if logmgr: + logmgr.close() + elif use_profiling: + print(actx.tabulate_profiling_data()) + + finish_tol = 1e-16 + assert np.abs(current_t - t_final) < finish_tol + + +if __name__ == "__main__": + import argparse + casename = "taylor-green-vortex" + parser = argparse.ArgumentParser(description=f"MIRGE-Com Example: {casename}") + parser.add_argument("--order", default=3, type=int, + help="specify the polynomial order of the DG method") + 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("--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("--numpy", action="store_true", + help="use numpy-based eager actx.") + parser.add_argument("--profiling", action="store_true", + help="turn on detailed performance profiling") + parser.add_argument("--log", action="store_true", default=True, + help="turn on logging") + parser.add_argument("--restart_file", help="root name of restart file") + parser.add_argument("--casename", help="casename to use for i/o") + args = parser.parse_args() + + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + + from mirgecom.array_context import get_reasonable_array_context_class + actx_class = get_reasonable_array_context_class( + lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) + + logging.basicConfig(format="%(message)s", level=logging.INFO) + if args.casename: + casename = args.casename + rst_filename = None + if args.restart_file: + rst_filename = args.restart_file + + main(actx_class, order=args.order, t_final=args.tfinal, + resolution=args.resolution, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, + casename=casename, rst_filename=rst_filename) + +# vim: foldmethod=marker diff --git a/examples/test_lazy_accuracy.sh b/examples/test_lazy_accuracy.sh index 5ba5785d5..25cc4c274 100755 --- a/examples/test_lazy_accuracy.sh +++ b/examples/test_lazy_accuracy.sh @@ -4,6 +4,8 @@ set -x set -e set -o pipefail +export POCL_DEBUG=cuda + test_list="vortex-mpi.py pulse-mpi.py" for file in ${test_list} do diff --git a/examples/thermally-coupled-mpi.py b/examples/thermally-coupled-mpi.py index 024fae806..8773dcb85 100644 --- a/examples/thermally-coupled-mpi.py +++ b/examples/thermally-coupled-mpi.py @@ -83,10 +83,8 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(use_logmgr=True, - use_overintegration=False, - use_leap=False, casename=None, - rst_filename=None, actx_class=None): +def main(actx_class, use_esdg=False, use_overintegration=False, + use_leap=False, casename=None, rst_filename=None): """Drive the example.""" if casename is None: casename = "mirgecom" @@ -99,7 +97,7 @@ def main(use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -175,7 +173,8 @@ def partition_generator_func(mesh, tag_to_elements, num_ranks): order = 3 dcoll = create_discretization_collection( - actx, volume_to_local_mesh, order=order) + actx, volume_to_local_mesh, order=order, + quadrature_order=order+2) dd_vol_fluid = DOFDesc(VolumeDomainTag("Fluid"), DISCR_TAG_BASE) dd_vol_wall = DOFDesc(VolumeDomainTag("Wall"), DISCR_TAG_BASE) @@ -504,8 +503,9 @@ def my_rhs(t, state, return_gradients=False): fluid_boundaries, wall_boundaries, fluid_state, wall_kappa, wall_temperature, time=t, - return_gradients=return_gradients, - quadrature_tag=quadrature_tag) + quadrature_tag=quadrature_tag, + use_esdg=use_esdg, + return_gradients=return_gradients) if return_gradients: ( @@ -572,8 +572,8 @@ def my_rhs_and_gradients(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") + parser.add_argument("--esdg", action="store_true", + help="use entropy-stable operator") parser.add_argument("--leap", action="store_true", help="use leap timestepper") parser.add_argument("--numpy", action="store_true", @@ -582,6 +582,14 @@ def my_rhs_and_gradients(t, state): parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -593,9 +601,9 @@ def my_rhs_and_gradients(t, state): if args.restart_file: rst_filename = args.restart_file - main(use_logmgr=args.log, use_overintegration=args.overintegration, + main(actx_class, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, use_leap=args.leap, - casename=casename, rst_filename=rst_filename, actx_class=actx_class, - ) + casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker diff --git a/examples/vortex-mpi.py b/examples/vortex-mpi.py index c2f86b9e9..c72b266e4 100644 --- a/examples/vortex-mpi.py +++ b/examples/vortex-mpi.py @@ -67,7 +67,7 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point -def main(actx_class, use_logmgr=True, +def main(actx_class, use_overintegration=False, use_esdg=False, use_leap=False, casename=None, rst_filename=None): """Drive the example.""" if casename is None: @@ -81,7 +81,7 @@ def main(actx_class, use_logmgr=True, from mirgecom.simutil import global_reduce as _global_reduce global_reduce = partial(_global_reduce, comm=comm) - logmgr = initialize_logmgr(use_logmgr, + logmgr = initialize_logmgr(True, filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) from mirgecom.array_context import initialize_actx, actx_class_is_profiling @@ -139,6 +139,12 @@ def main(actx_class, use_logmgr=True, dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) + from grudge.dof_desc import DISCR_TAG_QUAD + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None + vis_timer = None if logmgr: @@ -344,7 +350,8 @@ 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, - boundaries=boundaries, gas_model=gas_model) + boundaries=boundaries, gas_model=gas_model, + quadrature_tag=quadrature_tag, use_esdg=use_esdg) current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, t_final, constant_cfl) @@ -384,16 +391,26 @@ def my_rhs(t, state): help="switch to a lazy computation mode") parser.add_argument("--profiling", action="store_true", help="turn on detailed performance profiling") - parser.add_argument("--log", action="store_true", default=True, - help="turn on logging") parser.add_argument("--leap", action="store_true", help="use leap timestepper") + 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("--numpy", action="store_true", help="use numpy-based eager actx.") parser.add_argument("--restart_file", help="root name of restart file") parser.add_argument("--casename", help="casename to use for i/o") args = parser.parse_args() + from warnings import warn + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + from mirgecom.array_context import get_reasonable_array_context_class actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) @@ -405,7 +422,8 @@ def my_rhs(t, state): if args.restart_file: rst_filename = args.restart_file - main(actx_class, use_logmgr=args.log, use_leap=args.leap, - casename=casename, rst_filename=rst_filename) + main(actx_class, use_leap=args.leap, + casename=casename, rst_filename=rst_filename, + use_overintegration=args.overintegration or args.esdg, use_esdg=args.esdg) # vim: foldmethod=marker diff --git a/mirgecom/euler.py b/mirgecom/euler.py index ba2528ac7..c314a1260 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -63,21 +63,254 @@ ) from mirgecom.gas_model import make_operator_fluid_states -from mirgecom.inviscid import ( +from mirgecom.inviscid import ( # noqa inviscid_flux, inviscid_facial_flux_rusanov, - inviscid_flux_on_element_boundary + inviscid_flux_on_element_boundary, + entropy_stable_inviscid_facial_flux_rusanov, + entropy_stable_inviscid_facial_flux, + entropy_conserving_flux_chandrashekar, + entropy_conserving_flux_renac ) from mirgecom.operators import div_operator from mirgecom.utils import normalize_boundaries +from arraycontext import map_array_container +from mirgecom.gas_model import ( + project_fluid_state, + make_fluid_state_trace_pairs, + make_entropy_projected_fluid_state, + conservative_to_entropy_vars, + entropy_to_conservative_vars +) + +from meshmode.dof_array import DOFArray + +from functools import partial + +from grudge.trace_pair import ( + TracePair, + interior_trace_pairs, + tracepair_with_discr_tag +) +from grudge.projection import volume_quadrature_project +from grudge.flux_differencing import volume_flux_differencing + +import grudge.op as op + + +class _ESFluidCVTag(): + pass + + +class _ESFluidTemperatureTag(): + pass + + +def entropy_stable_euler_operator( + dcoll, gas_model, state, boundaries, time=0.0, + inviscid_numerical_flux_func=None, + entropy_conserving_flux_func=None, + operator_states_quad=None, + dd=DD_VOLUME_ALL, quadrature_tag=None, comm_tag=None, + limiter_func=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 + :class:`~grudge.dof_desc.BoundaryDomainTag` + + 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. + """ + boundaries = normalize_boundaries(boundaries) + + dd_vol = dd + dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) + dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) + + # 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_base = gas_model.eos.gamma(state.cv, state.temperature) + + # Interpolate state to vol quad grid + if operator_states_quad is not None: + state_quad = operator_states_quad[0] + else: + if state.is_mixture and limiter_func is None: + warn("Mixtures often require species limiting, and a non-limited " + "state is being created for this operator. For mixtures, " + "one should pass the operator_states_quad argument with " + "limited states, or least pass a limiter_func to this operator.") + state_quad = project_fluid_state( + dcoll, dd_vol, dd_vol_quad, state, gas_model, limiter_func=limiter_func, + entropy_stable=True) + + gamma_quad = gas_model.eos.gamma(state_quad.cv, state_quad.temperature) + + # Compute the projected (nodal) entropy variables + entropy_vars = volume_quadrature_project( + dcoll, dd_vol_quad, + # Map to entropy variables + conservative_to_entropy_vars(gamma_quad, state_quad)) + + modified_conserved_fluid_state = \ + make_entropy_projected_fluid_state(dcoll, dd_vol_quad, dd_allfaces_quad, + state, entropy_vars, gamma_base, + 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(dcoll.discr_from_dd(dd_vol).groups, ary))) + + if entropy_conserving_flux_func is None: + entropy_conserving_flux_func = \ + (entropy_conserving_flux_renac if state.is_mixture + else entropy_conserving_flux_chandrashekar) + flux_func = "renac" if state.is_mixture else "chandrashekar" + warn("No entropy_conserving_flux_func was given for ESDG. " + f"Setting EC flux to entropy_conserving_flux_{flux_func}.") + + flux_matrices = entropy_conserving_flux_func( + gas_model, + _reshape((1, -1), modified_conserved_fluid_state), + _reshape((-1, 1), modified_conserved_fluid_state)) + + # Compute volume derivatives using flux differencing + inviscid_vol_term = \ + -volume_flux_differencing(dcoll, dd_vol_quad, dd_allfaces_quad, + flux_matrices) + + # transfer trace pairs to quad grid, update pair dd + interp_to_surf_quad = partial(tracepair_with_discr_tag, dcoll, quadrature_tag) + + 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(dcoll, state.temperature, + volume_dd=dd_vol, + comm_tag=(_ESFluidTemperatureTag, + comm_tag)) + ] + + def _interp_to_surf_modified_conservedvars(gamma, ev_pair): + # Takes a trace pair containing the projected entropy variables + # and converts them into conserved variables on the quadrature grid. + local_dd = ev_pair.dd + local_dd_quad = local_dd.with_discr_tag(quadrature_tag) + + # Interpolate entropy variables to the surface quadrature grid + ev_pair_surf = op.project(dcoll, local_dd, local_dd_quad, ev_pair) + + if isinstance(gamma, DOFArray): + gamma = op.project(dcoll, dd_vol, local_dd_quad, gamma) + + return TracePair( + local_dd_quad, + # Convert interior and exterior states to conserved variables + interior=entropy_to_conservative_vars(gamma, ev_pair_surf.int), + exterior=entropy_to_conservative_vars(gamma, ev_pair_surf.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_base, tpair) + for tpair in interior_trace_pairs(dcoll, entropy_vars, volume_dd=dd_vol, + comm_tag=(_ESFluidCVTag, comm_tag))] + + boundary_states = { + # TODO: Use modified conserved vars as the input state? + # Would need to make an "entropy-projection" variant + # of *project_fluid_state* + bdtag: project_fluid_state( + dcoll, dd_vol, + # Make sure we get the state on the quadrature grid + # restricted to the tag *btag* + dd_vol_quad.with_domain_tag(bdtag), + state, gas_model, entropy_stable=True) for bdtag 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) + + if inviscid_numerical_flux_func is None: + inviscid_numerical_flux_func = \ + partial(entropy_stable_inviscid_facial_flux_rusanov, + entropy_conserving_flux_func=entropy_conserving_flux_func) + warn("No inviscid_numerical_flux_func was given for ESDG. " + "Automatically setting facial flux to entropy-stable Rusanov " + "(entropy_stable_inviscid_facial_flux_rusanov).") + elif inviscid_numerical_flux_func not in \ + [entropy_stable_inviscid_facial_flux_rusanov, + entropy_stable_inviscid_facial_flux]: + warn("Unrecognized inviscid_numerical_flux_func for ESDG. Proceed only " + "if you know what you are doing. An ESDG-compatible facial flux " + "function *must* be used with ESDG. Valid built-in choices are:\n" + "* entropy_stable_inviscid_facial_flux_rusanov, -or-\n" + "* entropy_stable_inviscid_facial_flux\n") + + # Compute interface contributions + inviscid_flux_bnd = inviscid_flux_on_element_boundary( + dcoll, gas_model, boundaries, interior_states, + boundary_states, quadrature_tag=quadrature_tag, + numerical_flux_func=inviscid_numerical_flux_func, time=time, + dd=dd_vol) + + return op.inverse_mass( + dcoll, + dd_vol, + inviscid_vol_term - op.face_mass(dcoll, dd_allfaces_quad, + inviscid_flux_bnd) + ) def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, - inviscid_numerical_flux_func=inviscid_facial_flux_rusanov, + inviscid_numerical_flux_func=None, quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, - comm_tag=None, - operator_states_quad=None): + comm_tag=None, use_esdg=False, operator_states_quad=None, + entropy_conserving_flux_func=None, limiter_func=None): r"""Compute RHS of the Euler flow equations. Returns @@ -138,14 +371,28 @@ def euler_operator(dcoll, state, gas_model, boundaries, time=0.0, dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) if operator_states_quad is None: - if state.is_mixture: + if state.is_mixture and limiter_func is None: warn("Mixtures often require species limiting, and a non-limited " "state is being created for this operator. For mixtures, " "one should pass the operator_states_quad argument with " - "limited states.") + "limited states or pass a limiter_func to this operator.") operator_states_quad = make_operator_fluid_states( dcoll, state, gas_model, boundaries, quadrature_tag, - dd=dd_vol, comm_tag=comm_tag) + dd=dd_vol, comm_tag=comm_tag, limiter_func=limiter_func, + entropy_stable=use_esdg) + + if use_esdg: + return entropy_stable_euler_operator( + dcoll, gas_model=gas_model, state=state, boundaries=boundaries, + time=time, operator_states_quad=operator_states_quad, dd=dd, + inviscid_numerical_flux_func=inviscid_numerical_flux_func, + entropy_conserving_flux_func=entropy_conserving_flux_func, + quadrature_tag=quadrature_tag, comm_tag=comm_tag) + + if inviscid_numerical_flux_func is None: + warn("inviscid_numerical_flux_func unspecified, defaulting to " + "inviscid_facial_flux_rusanov.") + inviscid_numerical_flux_func = inviscid_facial_flux_rusanov volume_state_quad, interior_state_pairs_quad, domain_boundary_states_quad = \ operator_states_quad diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py index 32f3e0e74..0de7d1501 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -19,6 +19,9 @@ .. autofunction:: project_fluid_state .. autofunction:: make_fluid_state_trace_pairs .. autofunction:: make_operator_fluid_states +.. autofunction:: make_entropy_projected_fluid_state +.. autofunction:: conservative_to_entropy_vars +.. autofunction:: entropy_to_conservative_vars """ __copyright__ = """ @@ -45,7 +48,7 @@ THE SOFTWARE. """ - +import numpy as np # noqa from functools import partial from dataclasses import dataclass from typing import Optional @@ -434,7 +437,8 @@ def make_fluid_state(cv, gas_model, return ViscousFluidState(cv=cv, dv=dv, tv=tv) -def project_fluid_state(dcoll, src, tgt, state, gas_model, limiter_func=None): +def project_fluid_state(dcoll, src, tgt, state, gas_model, limiter_func=None, + entropy_stable=False): """Project a fluid state onto a boundary consistent with the gas model. If required by the gas model, (e.g. gas is a mixture), this routine will @@ -483,6 +487,14 @@ def project_fluid_state(dcoll, src, tgt, state, gas_model, limiter_func=None): if state.is_mixture: temperature_seed = op.project(dcoll, src, tgt, state.dv.temperature) + if entropy_stable: + temp_state = make_fluid_state(cv=cv_sd, gas_model=gas_model, + temperature_seed=temperature_seed, + limiter_func=limiter_func, limiter_dd=tgt) + gamma = gas_model.eos.gamma(temp_state.cv, temp_state.temperature) + ev_sd = conservative_to_entropy_vars(gamma, temp_state) + cv_sd = entropy_to_conservative_vars(gamma, ev_sd) + smoothness_mu = None if state.dv.smoothness_mu is not None: smoothness_mu = op.project(dcoll, src, tgt, state.dv.smoothness_mu) @@ -622,7 +634,7 @@ class _WallDensityTag: def make_operator_fluid_states( dcoll, volume_state, gas_model, boundaries, quadrature_tag=DISCR_TAG_BASE, - dd=DD_VOLUME_ALL, comm_tag=None, limiter_func=None): + dd=DD_VOLUME_ALL, comm_tag=None, limiter_func=None, entropy_stable=False): """Prepare gas model-consistent fluid states for use in fluid operators. This routine prepares a model-consistent fluid state for each of the volume and @@ -695,9 +707,10 @@ def make_operator_fluid_states( interp_to_surf_quad = partial(tracepair_with_discr_tag, dcoll, quadrature_tag) domain_boundary_states_quad = { - bdtag: project_fluid_state(dcoll, dd_vol, - dd_vol_quad.with_domain_tag(bdtag), - volume_state, gas_model, limiter_func=limiter_func) + bdtag: project_fluid_state( + dcoll, dd_vol, dd_vol_quad.with_domain_tag(bdtag), + volume_state, gas_model, limiter_func=limiter_func, + entropy_stable=entropy_stable) for bdtag in boundaries } @@ -769,9 +782,9 @@ def make_operator_fluid_states( # Interpolate the fluid state to the volume quadrature grid # (this includes the conserved and dependent quantities) - volume_state_quad = project_fluid_state(dcoll, dd_vol, dd_vol_quad, - volume_state, gas_model, - limiter_func=limiter_func) + volume_state_quad = project_fluid_state( + dcoll, dd_vol, dd_vol_quad, volume_state, gas_model, + limiter_func=limiter_func, entropy_stable=entropy_stable) return \ volume_state_quad, interior_boundary_states_quad, domain_boundary_states_quad @@ -861,3 +874,113 @@ def replace_fluid_state( material_densities=material_densities, limiter_func=limiter_func, limiter_dd=limiter_dd) + + +def make_entropy_projected_fluid_state( + discr, dd_vol, dd_faces, state, entropy_vars, gamma, gas_model): + """Projects the entropy vars to target manifold, computes the CV from that.""" + 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, species densities) + into entropy variables, per eqn. 4.5 of [Chandrashekar_2013]_. + + 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 + y_species = state.species_mass_fractions + + u_square = np.dot(u, u) + s = actx.np.log(p) - gamma*actx.np.log(rho) + rho_p = rho / p + ev_mass = ((gamma - s) / (gamma - 1)) - 0.5 * rho_p * u_square + ev_spec = y_species / (gamma - 1) + + return make_conserved(dim, mass=ev_mass, energy=-rho_p, momentum=rho_p * u, + species_mass=ev_spec) + + +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, species density) per + eqn. 4.5 of [Chandrashekar_2013]_. + + 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 = np.dot(v234, v234) + s = gamma - v1 + v_square/(2*v5) + iota = ((gamma - 1) / (-v5)**gamma)**(inv_gamma_minus_one) + rho_iota = iota * actx.np.exp(-s * inv_gamma_minus_one) + mass = -rho_iota * v5 + spec_mass = mass * v6ns + + return make_conserved( + dim, + mass=mass, + energy=rho_iota * (1 - v_square/(2*v5)), + momentum=rho_iota * v234, + species_mass=spec_mass + ) diff --git a/mirgecom/initializers.py b/mirgecom/initializers.py index 040237b54..7eaaf4d43 100644 --- a/mirgecom/initializers.py +++ b/mirgecom/initializers.py @@ -15,6 +15,7 @@ .. autoclass:: ShearFlow .. autoclass:: PlanarDiscontinuity .. autoclass:: MulticomponentTrig +.. autoclass:: InviscidTaylorGreenVortex State Initializers ^^^^^^^^^^^^^^^^^^ @@ -261,7 +262,7 @@ class SodShock1D: """ def __init__( - self, *, dim=2, xdir=0, x0=0.5, rhol=1.0, + self, *, dim=1, xdir=0, x0=0.5, rhol=1.0, rhor=0.125, pleft=1.0, pright=0.1, ): """Initialize shock parameters. @@ -284,8 +285,8 @@ def __init__( self._x0 = x0 self._rhol = rhol self._rhor = rhor - self._energyl = pleft - self._energyr = pright + self._pl = pleft + self._pr = pright self._dim = dim self._xdir = xdir if self._xdir >= self._dim: @@ -304,30 +305,29 @@ def __call__(self, x_vec, *, eos=None, **kwargs): """ if eos is None: eos = IdealSingleGas() + gm1 = eos.gamma() - 1.0 gmn1 = 1.0 / gm1 - x_rel = x_vec[self._xdir] - actx = x_rel.array_context + x = x_vec[self._xdir] + actx = x.array_context - zeros = 0*x_rel + zeros = actx.np.zeros_like(x) rhor = zeros + self._rhor rhol = zeros + self._rhol x0 = zeros + self._x0 - energyl = zeros + gmn1 * self._energyl - energyr = zeros + gmn1 * self._energyr - yesno = actx.np.greater(x_rel, x0) - mass = actx.np.where(yesno, rhor, rhol) - energy = actx.np.where(yesno, energyr, energyl) - mom = make_obj_array( - [ - 0*x_rel - for i in range(self._dim) - ] - ) + energyl = zeros + gmn1 * self._pl + energyr = zeros + gmn1 * self._pr + + sigma = 1e-13 + weight = 0.5 * (1.0 - actx.np.tanh(1.0/sigma * (x - x0))) + + mass = rhor + (rhol - rhor)*weight + energy = energyr + (energyl - energyr)*weight + momentum = make_obj_array([1.*zeros for _ in range(self._dim)]) return make_conserved(dim=self._dim, mass=mass, energy=energy, - momentum=mom) + momentum=momentum) class DoubleMachReflection: @@ -818,7 +818,8 @@ def __init__( spec_centers=None, spec_omegas=None, spec_diffusivities=None, - wave_vector=None + wave_vector=None, + trig_function=None ): r"""Initialize MulticomponentLump parameters. @@ -835,6 +836,10 @@ def __init__( velocity: numpy.ndarray fixed flow velocity used for exact solution at t != 0, shape ``(dim,)`` + wave_vector: numpy.ndarray + optional fixed vector indicating normal direction of wave + trig_function + callable trig function """ if center is None: center = np.zeros(shape=(dim,)) @@ -843,7 +848,7 @@ def __init__( if center.shape != (dim,) or velocity.shape != (dim,): raise ValueError(f"Expected {dim}-dimensional vector inputs.") if spec_y0s is None: - spec_y0s = np.ones(shape=(nspecies,)) + spec_y0s = 2.0*np.ones(shape=(nspecies,)) if spec_centers is None: spec_centers = make_obj_array([np.zeros(shape=dim,) for i in range(nspecies)]) @@ -852,6 +857,7 @@ def __init__( if spec_amplitudes is None: spec_amplitudes = np.ones(shape=(nspecies,)) + if spec_diffusivities is None: spec_diffusivities = np.ones(shape=(nspecies,)) @@ -859,6 +865,10 @@ def __init__( wave_vector = np.zeros(shape=(dim,)) wave_vector[0] = 1 + import mirgecom.math as mm + if trig_function is None: + trig_function = mm.sin + if len(spec_y0s) != nspecies or\ len(spec_amplitudes) != nspecies or\ len(spec_centers) != nspecies: @@ -881,6 +891,7 @@ def __init__( self._spec_omegas = spec_omegas self._d = spec_diffusivities self._wave_vector = wave_vector + self._trig_func = trig_function def __call__(self, x_vec, *, eos=None, time=0, **kwargs): """ @@ -910,15 +921,15 @@ def __call__(self, x_vec, *, eos=None, time=0, **kwargs): energy = ((self._p0 / (self._gamma - 1.0)) + 0.5*mass*np.dot(self._velocity, self._velocity)) - import mirgecom.math as mm vel_t = t * self._velocity + import mirgecom.math as mm spec_mass = np.empty((self._nspecies,), dtype=object) for i in range(self._nspecies): spec_x = x_vec - self._spec_centers[i] wave_r = spec_x - vel_t wave_x = np.dot(wave_r, self._wave_vector) expterm = mm.exp(-t*self._d[i]*self._spec_omegas[i]**2) - trigterm = mm.sin(self._spec_omegas[i]*wave_x) + trigterm = self._trig_func(self._spec_omegas[i]*wave_x) spec_y = self._spec_y0s[i] + self._spec_amps[i]*expterm*trigterm spec_mass[i] = mass * spec_y @@ -1552,3 +1563,74 @@ def __call__(self, x, **kwargs): return make_conserved(dim=self._dim, mass=density, momentum=mom, energy=total_energy) + + +class InviscidTaylorGreenVortex: + """Initialize Taylor-Green Vortex.""" + + def __init__( + self, *, dim=3, mach_number=0.05, domain_lengthscale=1, v0=1, p0=1, + viscosity=1e-5 + ): + """Initialize vortex parameters.""" + self._mach_number = mach_number + self._domain_lengthscale = domain_lengthscale + self._v0 = v0 + self._p0 = p0 + self._mu = viscosity + self._dim = dim + + def __call__(self, x_vec, *, eos=None, time=0, **kwargs): + """ + Create the 3D Taylor-Green initial profile at locations *x_vec*. + + Parameters + ---------- + x_vec: numpy.ndarray + Nodal coordinates + eos: :class:`mirgecom.eos.IdealSingleGas` + Equation of state class with method to supply gas *gamma*. + """ + if eos is None: + eos = IdealSingleGas() + + length = self._domain_lengthscale + gamma = eos.gamma() + v0 = self._v0 + p0 = self._p0 + rho0 = gamma * self._mach_number ** 2 + dim = len(x_vec) + x = x_vec[0] + y = x_vec[1] + actx = x_vec[0].array_context + zeros = actx.zeros_like(x) + ones = 1 + zeros + nu = self._mu/rho0 + ft = actx.np.exp(-2*nu*time) + + if dim == 3: + z = x_vec[0] + + p = p0 + rho0 * (v0 ** 2) / 16 * ( + actx.np.cos(2*x / length + actx.np.cos(2*y / length)) + ) * actx.np.cos(2*z / length + 2) + u = ( + v0 * actx.np.sin(x / length) * actx.np.cos(y / length) + ) * actx.np.cos(z / length) + v = ( + -v0 * actx.np.cos(x / length) * actx.np.sin(y / length) + ) * actx.np.cos(z / length) + w = zeros + velocity = make_obj_array([u, v, w]) + else: + u = actx.np.sin(x)*actx.np.cos(y)*ft + v = -actx.np.cos(x)*actx.np.sin(y)*ft + p = rho0/4.0 * (actx.np.cos(2*x) + actx.np.sin(2*y)) * ft * ft + velocity = make_obj_array([u, v]) + + momentum = rho0 * velocity + energy = p / (gamma - 1) + rho0 / 2 * np.dot(velocity, velocity) + rho = rho0 * ones + + return make_conserved(dim=self._dim, mass=rho, + energy=energy, momentum=momentum) diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index c245ee635..f1c713a6d 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -7,6 +7,10 @@ .. autofunction:: inviscid_facial_flux_rusanov .. autofunction:: inviscid_facial_flux_hll .. autofunction:: inviscid_flux_on_element_boundary +.. autofunction:: entropy_conserving_flux_chandrashekar +.. autofunction:: entropy_conserving_flux_renac +.. autofunction:: entropy_stable_inviscid_facial_flux +.. autofunction:: entropy_stable_inviscid_facial_flux_rusanov Inviscid Time Step Computation ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -47,9 +51,15 @@ DISCR_TAG_BASE, ) import grudge.op as op -from mirgecom.fluid import make_conserved +from mirgecom.fluid import ( + make_conserved, + ConservedVars +) from mirgecom.utils import normalize_boundaries +from arraycontext import outer +from meshmode.dof_array import DOFArray + def inviscid_flux(state): r"""Compute the inviscid flux vectors from fluid conserved vars *cv*. @@ -388,3 +398,266 @@ def get_inviscid_cfl(dcoll, state, dt): The CFL at each node. """ return dt / get_inviscid_timestep(dcoll, state=state) + + +def entropy_conserving_flux_chandrashekar(gas_model, state_ll, state_rr): + """Compute the entropy conservative fluxes from states *state_ll* and *state_rr*. + + This routine implements the two-point volume flux based on the entropy + conserving and kinetic energy preserving two-point flux in + equations (4.12 - 4.14) of [Chandrashekar_2013]_. + + 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 + y_ll = state_ll.species_mass_fractions + + rho_rr = state_rr.mass_density + u_rr = state_rr.velocity + p_rr = state_rr.pressure + y_rr = state_rr.species_mass_fractions + + beta_ll = 0.5 * rho_ll / p_ll + beta_rr = 0.5 * rho_rr / p_rr + specific_kin_ll = 0.5 * np.dot(u_ll, u_ll) + specific_kin_rr = 0.5 * np.dot(u_rr, u_rr) + + rho_avg = 0.5 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) + y_mean = 0.5 * (y_ll + y_rr) + rho_species_mean = rho_mean * y_mean + + 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_conserving_flux_renac(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 + equation (24) of [Renac_2021]_. + + 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 + t_ll = state_ll.temperature + t_rr = state_rr.temperature + p_ll = state_ll.pressure + p_rr = state_rr.pressure + gamma_ll = gas_model.eos.gamma(state_ll.cv, state_ll.temperature) + gamma_rr = gas_model.eos.gamma(state_rr.cv, state_rr.temperature) + theta_ll = 1.0/t_ll + theta_rr = 1.0/t_rr + t_avg = 0.5*(t_ll + t_rr) + + pot_ll = p_ll * theta_ll + pot_rr = p_rr * theta_rr + pot_avg = 0.5*(pot_ll + pot_rr) + + 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) + ) + + theta_mean = ln_mean(theta_ll, theta_rr) + t_mean = 1.0/theta_mean + pec_avg = pot_avg * t_avg + p_mean = ln_mean(p_ll, p_rr) + + # Primitive variables for left and right states + rho_ll = state_ll.mass_density + u_ll = state_ll.velocity + p_ll = state_ll.pressure + y_ll = state_ll.species_mass_fractions + + rho_rr = state_rr.mass_density + u_rr = state_rr.velocity + p_rr = state_rr.pressure + y_rr = state_rr.species_mass_fractions + + kin_comb = 0.5 * np.dot(u_ll, u_rr) + + rho_mean = ln_mean(rho_ll, rho_rr) + y_avg = 0.5 * (y_ll + y_rr) + species_mass_mean = rho_mean * y_avg + + u_avg = 0.5 * (u_ll + u_rr) + + mass_flux = rho_mean * u_avg + momentum_flux = outer(mass_flux, u_avg) + np.eye(dim) * pec_avg + + gamma_avg = 0.5 * (gamma_ll + gamma_rr) + ener_es = p_mean / (gamma_avg - 1) + 0.5 * rho_mean * np.dot(u_avg, u_avg) + cv_es = ConservedVars(mass=rho_mean, momentum=rho_mean*u_avg, + species_mass=species_mass_mean, energy=ener_es) + heat_cap_cv_es_mix = gas_model.eos.heat_capacity_cv(cv_es, t_mean) + ener_term = (heat_cap_cv_es_mix * t_mean + kin_comb) * mass_flux + energy_flux = ener_term + pec_avg * u_avg + + species_mass_flux = species_mass_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_facial_flux(state_pair, gas_model, normal, + entropy_conserving_flux_func=None, + alpha=None): + r"""Return the entropy stable inviscid numerical flux across a face. + + 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 optionally imposed by penalizing the "jump" of the state across + interfaces with strength *alpha*. + + 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 + + gas_model: :class:`~mirgecom.gas_model.GasModel` + + Physical gas model including equation of state, transport, + and kinetic properties as required by fluid state + + normal: numpy.ndarray + + The element interface normals + + entropy_conserving_flux_func: + + Callable function returning the entropy-conserving flux function to + use. If unspecified, an appropriate flux function will be chosen + depending on the type of fluid state (e.g. mixture vs. single gas). + + alpha: + + Strength of the penalization term. This can be a fixed single scalar, + or a :class:`~meshmode.dof_array.DOFArray`. For example, a Rusanov flux can + be constructed by passing the max wavespeed as alpha. + + Returns + ------- + :class:`~mirgecom.fluid.ConservedVars` + A CV object containing the scalar numerical fluxes at the input faces. + """ + # Automatically choose the appropriate EC flux if none is given + if entropy_conserving_flux_func is None: + entropy_conserving_flux_func = \ + (entropy_conserving_flux_renac if state_pair.int.is_mixture + else entropy_conserving_flux_chandrashekar) + if state_pair.int.is_mixture: + from warnings import warn + warn("`entropy_conserving_flux_renac` is expensive to compile for " + "mixtures.") + + flux = entropy_conserving_flux_func(gas_model, + state_pair.int, + state_pair.ext) + + if alpha is not None: + flux = flux - 0.5*alpha*outer(state_pair.ext.cv - state_pair.int.cv, normal) + + return flux @ normal + + +def entropy_stable_inviscid_facial_flux_rusanov(state_pair, gas_model, normal, + entropy_conserving_flux_func=None, + **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`). Rusanov + dissipation is imposed by penalizing the "jump" of the state across + interfaces with the max wavespeed between the two (+/-) facial states. + + 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 + + gas_model: :class:`~mirgecom.gas_model.GasModel` + + Physical gas model including equation of state, transport, + and kinetic properties as required by fluid state + + normal: numpy.ndarray + + The element interface normals + + entropy_conserving_flux_func: + + Callable function returning the entropy-conserving flux function to + use. If unspecified, an appropriate flux function will be chosen + depending on the type of fluid state (e.g. mixture vs. single gas). + + Returns + ------- + :class:`~mirgecom.fluid.ConservedVars` + A CV object containing the scalar numerical fluxes at the input faces. + """ + # This calculates the local maximum eigenvalue of the flux Jacobian + # for a single component gas, i.e. the element-local max wavespeed |v| + c. + actx = state_pair.int.array_context + alpha = actx.np.maximum(state_pair.int.wavespeed, state_pair.ext.wavespeed) + + return entropy_stable_inviscid_facial_flux( + state_pair=state_pair, gas_model=gas_model, normal=normal, + entropy_conserving_flux_func=entropy_conserving_flux_func, + alpha=alpha) diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py index 171c41555..f8e411387 100644 --- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py +++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py @@ -91,7 +91,6 @@ IsothermalSlipWallBoundary, IsothermalWallBoundary) from mirgecom.flux import num_flux_central -from mirgecom.inviscid import inviscid_facial_flux_rusanov from mirgecom.viscous import viscous_facial_flux_harmonic from mirgecom.gas_model import ( replace_fluid_state, @@ -1375,9 +1374,8 @@ def coupled_ns_heat_operator( quadrature_tag=DISCR_TAG_BASE, limiter_func=None, fluid_gradient_numerical_flux_func=num_flux_central, - inviscid_numerical_flux_func=inviscid_facial_flux_rusanov, - viscous_numerical_flux_func=viscous_facial_flux_harmonic, - return_gradients=False): + return_gradients=False, + ns_operator=ns_operator): r""" Compute the RHS of the fluid and wall subdomains. @@ -1467,18 +1465,6 @@ def coupled_ns_heat_operator( the temperature gradient in the fluid subdomain. Defaults to :class:`~mirgecom.flux.num_flux_central`. - inviscid_numerical_flux_func: - - Callable function providing the face-normal flux to be used - for the divergence of the inviscid transport flux. This defaults to - :func:`~mirgecom.inviscid.inviscid_facial_flux_rusanov`. - - viscous_numerical_flux_func: - - Callable function providing the face-normal flux to be used - for the divergence of the viscous transport flux. This defaults to - :func:`~mirgecom.viscous.viscous_facial_flux_harmonic`. - limiter_func: Callable function to be passed to @@ -1570,11 +1556,7 @@ def coupled_ns_heat_operator( quadrature_tag=quadrature_tag) # Compute the subdomain NS/diffusion operators using the augmented boundaries - - my_ns_operator = partial(ns_operator, - inviscid_numerical_flux_func=inviscid_numerical_flux_func, - viscous_numerical_flux_func=viscous_numerical_flux_func) - ns_result = my_ns_operator( + ns_result = ns_operator( dcoll, gas_model, fluid_state, fluid_all_boundaries, time=time, quadrature_tag=quadrature_tag, dd=fluid_dd, return_gradients=return_gradients, @@ -1613,7 +1595,8 @@ def basic_coupled_ns_heat_operator( wall_penalty_amount=None, quadrature_tag=DISCR_TAG_BASE, limiter_func=None, - return_gradients=False): + return_gradients=False, + use_esdg=False): r""" Simple implementation of a thermally-coupled fluid/wall operator. @@ -1712,6 +1695,10 @@ def basic_coupled_ns_heat_operator( that filters or limits the produced fluid states. This is used to keep species mass fractions in physical and realizable states, for example. + use_esdg: bool + + If `True`, use the entropy-stable version of the Navier-Stokes operator. + Returns ------- @@ -1788,7 +1775,8 @@ def basic_coupled_ns_heat_operator( # Compute the subdomain NS/diffusion operators using the augmented boundaries my_ns_operator = partial(ns_operator, - viscous_numerical_flux_func=viscous_facial_flux_harmonic) + viscous_numerical_flux_func=viscous_facial_flux_harmonic, + use_esdg=use_esdg) ns_result = my_ns_operator( dcoll, gas_model, fluid_state, fluid_all_boundaries, time=time, quadrature_tag=quadrature_tag, dd=fluid_dd, diff --git a/mirgecom/navierstokes.py b/mirgecom/navierstokes.py index 64728b428..fcb095406 100644 --- a/mirgecom/navierstokes.py +++ b/mirgecom/navierstokes.py @@ -75,10 +75,14 @@ import grudge.op as op +from mirgecom.euler import ( + euler_operator, + entropy_stable_euler_operator +) from mirgecom.inviscid import ( inviscid_flux, inviscid_facial_flux_rusanov, - inviscid_flux_on_element_boundary + inviscid_flux_on_element_boundary, ) from mirgecom.viscous import ( viscous_flux, @@ -102,6 +106,14 @@ class _NSGradTemperatureTag: pass +class _ESFluidCVTag(): + pass + + +class _ESFluidTemperatureTag(): + pass + + def _gradient_flux_interior(dcoll, numerical_flux_func, tpair): """Compute interior face flux for gradient operator.""" from arraycontext import outer @@ -119,7 +131,8 @@ def grad_cv_operator( quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, comm_tag=None, # Added to avoid repeated computation # FIXME: See if there's a better way to do this - operator_states_quad=None): + operator_states_quad=None, limiter_func=None, + use_esdg=False): r"""Compute the gradient of the fluid conserved variables. Parameters @@ -176,14 +189,15 @@ def grad_cv_operator( dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) if operator_states_quad is None: - if state.is_mixture: + if state.is_mixture and limiter_func is None: warn("Mixtures often require species limiting, and a non-limited " "state is being created for this operator. For mixtures, " "one should pass the operator_states_quad argument with " - "limited states.") + "limited states or pass a limiter_func to this operator.") operator_states_quad = make_operator_fluid_states( dcoll, state, gas_model, boundaries, quadrature_tag, - dd=dd_vol, comm_tag=comm_tag) + dd=dd_vol, comm_tag=comm_tag, limiter_func=limiter_func, + entropy_stable=use_esdg) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -226,7 +240,7 @@ def grad_t_operator( quadrature_tag=DISCR_TAG_BASE, dd=DD_VOLUME_ALL, comm_tag=None, # Added to avoid repeated computation # FIXME: See if there's a better way to do this - operator_states_quad=None): + operator_states_quad=None, limiter_func=None, use_esdg=False): r"""Compute the gradient of the fluid temperature. Parameters @@ -282,14 +296,15 @@ def grad_t_operator( dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) if operator_states_quad is None: - if state.is_mixture: + if state.is_mixture and limiter_func is None: warn("Mixtures often require species limiting, and a non-limited " "state is being created for this operator. For mixtures, " "one should pass the operator_states_quad argument with " - "limited states.") + "limited states or pass a limiter_func to this operator.") operator_states_quad = make_operator_fluid_states( dcoll, state, gas_model, boundaries, quadrature_tag, - dd=dd_vol, comm_tag=comm_tag) + dd=dd_vol, comm_tag=comm_tag, limiter_func=limiter_func, + entropy_stable=use_esdg) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -331,15 +346,17 @@ def grad_t_operator( def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, + inviscid_fluid_operator=None, inviscid_numerical_flux_func=inviscid_facial_flux_rusanov, gradient_numerical_flux_func=num_flux_central, viscous_numerical_flux_func=viscous_facial_flux_central, return_gradients=False, quadrature_tag=DISCR_TAG_BASE, - dd=DD_VOLUME_ALL, comm_tag=None, + dd=DD_VOLUME_ALL, comm_tag=None, limiter_func=None, # Added to avoid repeated computation # FIXME: See if there's a better way to do this - operator_states_quad=None, - grad_cv=None, grad_t=None, inviscid_terms_on=True): + operator_states_quad=None, use_esdg=False, + grad_cv=None, grad_t=None, inviscid_terms_on=True, + entropy_conserving_flux_func=None): r"""Compute RHS of the Navier-Stokes equations. Parameters @@ -432,6 +449,11 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, if dd.discretization_tag != DISCR_TAG_BASE: raise ValueError("dd must belong to the base discretization") + if inviscid_fluid_operator is None and use_esdg: + inviscid_fluid_operator = entropy_stable_euler_operator + elif use_esdg and inviscid_fluid_operator == euler_operator: + raise RuntimeError("Standard Euler operator is incompatible with ESDG.") + dd_vol = dd dd_vol_quad = dd_vol.with_discr_tag(quadrature_tag) dd_allfaces_quad = dd_vol_quad.trace(FACE_RESTR_ALL) @@ -444,14 +466,14 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, # Note: these states will live on the quadrature domain if one is given, # otherwise they stay on the interpolatory/base domain. if operator_states_quad is None: - if state.is_mixture: + if state.is_mixture and limiter_func is None: warn("Mixtures often require species limiting, and a non-limited " "state is being created for this operator. For mixtures, " "one should pass the operator_states_quad argument with " - "limited states.") + "limited states or provide a limiter_func to this operator.") operator_states_quad = make_operator_fluid_states( dcoll, state, gas_model, boundaries, quadrature_tag, - dd=dd_vol, comm_tag=comm_tag) + limiter_func=limiter_func, dd=dd_vol, comm_tag=comm_tag) vol_state_quad, inter_elem_bnd_states_quad, domain_bnd_states_quad = \ operator_states_quad @@ -521,7 +543,10 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, dd=dd_vol) # Add corresponding inviscid parts if enabled - if inviscid_terms_on: + # Note that this is the default, and highest performing path. + # To get separate inviscid operator, set inviscid_fluid_operator explicitly + # or use ESDG. + if inviscid_terms_on and inviscid_fluid_operator is None: vol_term = vol_term - inviscid_flux(state=vol_state_quad) bnd_term = bnd_term - inviscid_flux_on_element_boundary( dcoll, gas_model, boundaries, inter_elem_bnd_states_quad, @@ -531,6 +556,14 @@ def ns_operator(dcoll, gas_model, state, boundaries, *, time=0.0, ns_rhs = div_operator(dcoll, dd_vol_quad, dd_allfaces_quad, vol_term, bnd_term) + # Call an external operator for the inviscid terms (Euler by default) + # ESDG *always* uses this branch + if inviscid_terms_on and inviscid_fluid_operator is not None: + ns_rhs = ns_rhs + inviscid_fluid_operator( + dcoll, state=state, gas_model=gas_model, boundaries=boundaries, + time=time, dd=dd, comm_tag=comm_tag, quadrature_tag=quadrature_tag, + operator_states_quad=operator_states_quad) + if return_gradients: return ns_rhs, grad_cv, grad_t return ns_rhs diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index a8da2609e..1c850be8b 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -38,6 +38,13 @@ .. autofunction:: compare_files_vtu .. autofunction:: compare_files_xdmf .. autofunction:: compare_files_hdf5 + +Exceptions +---------- + +.. autoclass:: SimulationConfigurationError +.. autoclass:: ApplicationOptionsError +.. autoclass:: SimulationRuntimeError """ __copyright__ = """ @@ -86,6 +93,24 @@ from mpi4py.MPI import Comm +class SimulationConfigurationError(RuntimeError): + """Simulation physics configuration or parameters error.""" + + pass + + +class SimulationRuntimeError(RuntimeError): + """General simulation runtime error.""" + + pass + + +class ApplicationOptionsError(RuntimeError): + """Application command-line options error.""" + + pass + + def get_number_of_tetrahedron_nodes(dim, order, include_faces=False): """Get number of nodes (modes) in *dim* Tetrahedron of *order*.""" # number of {nodes, modes} see e.g.: