diff --git a/MANIFEST.in b/MANIFEST.in index 5b36d2cae..9e8eebd39 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,2 +1,4 @@ include mirgecom/mechanisms/*.yaml include mirgecom/mechanisms/*.cti +include examples/*.geo +include examples/*.msh diff --git a/examples/doublemach_physical_av-mpi-lazy.py b/examples/doublemach_physical_av-mpi-lazy.py new file mode 120000 index 000000000..a974dc919 --- /dev/null +++ b/examples/doublemach_physical_av-mpi-lazy.py @@ -0,0 +1 @@ +doublemach_physical_av-mpi.py \ No newline at end of file diff --git a/examples/doublemach_physical_av-mpi.py b/examples/doublemach_physical_av-mpi.py new file mode 100644 index 000000000..d667da451 --- /dev/null +++ b/examples/doublemach_physical_av-mpi.py @@ -0,0 +1,750 @@ +"""Demonstrate double mach reflection.""" + +__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 +import pyopencl as cl +from functools import partial + +from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from grudge.dof_desc import BoundaryDomainTag +from grudge.shortcuts import make_visualizer + +from mirgecom.euler import euler_operator +from mirgecom.navierstokes import ns_operator +from mirgecom.artificial_viscosity import ( + av_laplacian_operator, + smoothness_indicator +) +from mirgecom.io import make_init_message +from mirgecom.mpi import mpi_entry_point +from mirgecom.integrators import rk4_step, euler_step +from grudge.shortcuts import compiled_lsrk45_step +from mirgecom.steppers import advance_state +from mirgecom.boundary import ( + PrescribedFluidBoundary, + PressureOutflowBoundary, + AdiabaticSlipBoundary +) +from mirgecom.initializers import DoubleMachReflection +from mirgecom.eos import IdealSingleGas +from mirgecom.transport import ( + SimpleTransport, + ArtificialViscosityTransportDiv +) +from mirgecom.simutil import get_sim_timestep +from mirgecom.utils import force_evaluation + +from logpyle import set_dt +from mirgecom.logging_quantities import ( + initialize_logmgr, + logmgr_add_cl_device_info, + logmgr_add_device_memory_usage +) + +logger = logging.getLogger(__name__) + + +class MyRuntimeError(RuntimeError): + """Simple exception to kill the simulation.""" + + pass + + +def get_doublemach_mesh(): + """Generate or import a grid using `gmsh`. + + Input required: + doubleMach.msh (read existing mesh) + + This routine will generate a new grid if it does + not find the grid file (doubleMach.msh). + """ + from meshmode.mesh.io import ( + read_gmsh, + generate_gmsh, + ScriptSource, + ) + import os + meshfile = "doubleMach.msh" + if not os.path.exists(meshfile): + mesh = generate_gmsh( + ScriptSource(""" + x0=1.0/6.0; + setsize=0.025; + Point(1) = {0, 0, 0, setsize}; + Point(2) = {x0,0, 0, setsize}; + Point(3) = {4, 0, 0, setsize}; + Point(4) = {4, 1, 0, setsize}; + Point(5) = {0, 1, 0, setsize}; + Line(1) = {1, 2}; + Line(2) = {2, 3}; + Line(5) = {3, 4}; + Line(6) = {4, 5}; + Line(7) = {5, 1}; + Line Loop(8) = {-5, -6, -7, -1, -2}; + Plane Surface(8) = {8}; + Physical Surface('domain') = {8}; + Physical Curve('flow') = {1, 6, 7}; + Physical Curve('wall') = {2}; + Physical Curve('out') = {5}; + """, "geo"), force_ambient_dim=2, dimensions=2, target_unit="M", + output_file_name=meshfile) + else: + mesh = read_gmsh(meshfile, force_ambient_dim=2) + + return mesh + + +@mpi_entry_point +def main(ctx_factory=cl.create_some_context, use_logmgr=True, + use_leap=False, use_profiling=False, use_overintegration=False, + casename=None, rst_filename=None, actx_class=None, lazy=False): + """Drive the example.""" + if actx_class is None: + raise RuntimeError("Array context class missing.") + + cl_ctx = ctx_factory() + + if casename is None: + casename = "mirgecom" + + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nparts = comm.Get_size() + + # logging and profiling + log_path = "log_data/" + logname = log_path + casename + ".sqlite" + + if rank == 0: + import os + log_dir = os.path.dirname(logname) + if log_dir and not os.path.exists(log_dir): + os.makedirs(log_dir) + comm.Barrier() + + logmgr = initialize_logmgr(use_logmgr, + filename=logname, mode="wo", mpi_comm=comm) + + if use_profiling: + queue = cl.CommandQueue( + cl_ctx, properties=cl.command_queue_properties.PROFILING_ENABLE) + else: + queue = cl.CommandQueue(cl_ctx) + + from mirgecom.simutil import get_reasonable_memory_pool + alloc = get_reasonable_memory_pool(cl_ctx, queue) + + if lazy: + actx = actx_class(comm, queue, mpi_base_tag=12000, allocator=alloc) + else: + actx = actx_class(comm, queue, allocator=alloc, force_device_scalars=True) + + # Timestepping control + current_step = 0 + timestepper = rk4_step + timestepper = euler_step + force_eval = True + t_final = 5.e-4 + current_cfl = 0.1 + current_dt = 2.5e-5 + current_t = 0 + constant_cfl = False + + def _compiled_stepper_wrapper(state, t, dt, rhs): + return compiled_lsrk45_step(actx, state, t, dt, rhs) + + timestepper = _compiled_stepper_wrapper + force_eval = False + + # default health status bounds + health_pres_min = 0.7 + health_pres_max = 19. + health_temp_min = 1e-4 + health_temp_max = 100 + + # Some i/o frequencies + nviz = 250 + nrestart = 1000 + nhealth = 1 + nstatus = 1 + + viz_path = "viz_data/" + vizname = viz_path + casename + 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 + gen_grid = partial(get_doublemach_mesh) + from mirgecom.simutil import generate_and_distribute_mesh + local_mesh, global_nelements = generate_and_distribute_mesh(comm, + gen_grid) + local_nelements = local_mesh.nelements + + from mirgecom.simutil import global_reduce as _global_reduce + global_reduce = partial(_global_reduce, comm=comm) + + from mirgecom.discretization import create_discretization_collection + order = 3 + 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 # noqa + + dim = 2 + if logmgr: + logmgr_add_cl_device_info(logmgr, queue) + logmgr_add_device_memory_usage(logmgr, queue) + + logmgr.add_watches([ + ("step.max", "step = {value}, "), + ("t_sim.max", "sim time: {value:1.6e} s\n"), + ("t_step.max", "------- step walltime: {value:6g} s, ") + ]) + + # which kind of artificial viscosity? + # 0 - none + # 1 - laplacian diffusion + # 2 - physical viscosity based, div(velocity) indicator + use_av = 2 + + # Solution setup and initialization + # {{{ Initialize simple transport model + kappa = 1.0 + s0 = np.log10(1.0e-4 / np.power(order, 4)) + alpha = 0.03 + kappa_t = 0. + sigma_v = 0. + # }}} + + # Shock strength + shock_location = 1.0/6.0 + shock_speed = 4.0 + shock_sigma = 0.01 + + initializer = DoubleMachReflection(shock_location=shock_location, + shock_speed=shock_speed, + shock_sigma=shock_sigma) + + from mirgecom.gas_model import GasModel, make_fluid_state + physical_transport = SimpleTransport( + viscosity=sigma_v, thermal_conductivity=kappa_t) + if use_av < 2: + transport_model = physical_transport + else: + transport_model = ArtificialViscosityTransportDiv( + physical_transport=physical_transport, + av_mu=1.0, av_prandtl=0.75) + + eos = IdealSingleGas() + gas_model = GasModel(eos=eos, transport=transport_model) + + def get_fluid_state(cv, smoothness_mu=None): + return make_fluid_state(cv=cv, gas_model=gas_model, + smoothness_mu=smoothness_mu) + + create_fluid_state = actx.compile(get_fluid_state) + + from grudge.dt_utils import characteristic_lengthscales + length_scales = characteristic_lengthscales(actx, dcoll) + + from mirgecom.navierstokes import grad_cv_operator + + # compiled wrapper for grad_cv_operator + def _grad_cv_operator(fluid_state, time): + return grad_cv_operator(dcoll=dcoll, gas_model=gas_model, + boundaries=boundaries, + state=fluid_state, + time=time, + quadrature_tag=quadrature_tag) + + grad_cv_operator_compiled = actx.compile(_grad_cv_operator) # noqa + + def compute_smoothness(cv, grad_cv): + + actx = cv.array_context + from mirgecom.fluid import velocity_gradient + div_v = np.trace(velocity_gradient(cv, grad_cv)) + + # kappa_h = 1.5 + kappa_h = 5 + gamma = gas_model.eos.gamma(cv) + r = gas_model.eos.gas_const(cv) + static_temp = 0.015 + c_star = actx.np.sqrt(gamma*r*(2/(gamma+1)*static_temp)) + indicator = -kappa_h*length_scales*div_v/c_star + + # steepness of the smoothed function + alpha = 100 + # cutoff, smoothness below this value is ignored + beta = 0.01 + smoothness = actx.np.log( + 1 + actx.np.exp(alpha*(indicator - beta)))/alpha + return smoothness*kappa_h*length_scales + + compute_smoothness_compiled = actx.compile(compute_smoothness) # noqa + + 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 = initializer(nodes) + + smoothness = None + no_smoothness = None + if use_av > 0: + smoothness = smoothness_indicator(dcoll, current_cv.mass, + kappa=kappa, s0=s0) + no_smoothness = 0.*smoothness + + current_state = make_fluid_state(cv=current_cv, gas_model=gas_model, + smoothness_mu=smoothness) + force_evaluation(actx, current_state) + + def _boundary_state(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(cv=initializer(x_vec=nodes, eos=gas_model.eos, + **kwargs), + gas_model=gas_model, + smoothness_mu=state_minus.dv.smoothness_mu) + + flow_boundary = PrescribedFluidBoundary( + boundary_state_func=_boundary_state) + + boundaries = { + BoundaryDomainTag("flow"): flow_boundary, + BoundaryDomainTag("wall"): AdiabaticSlipBoundary(), + BoundaryDomainTag("out"): PressureOutflowBoundary(boundary_pressure=1.0), + } + + visualizer = make_visualizer(dcoll, order) + + initname = initializer.__class__.__name__ + 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 vol_min(x): + from grudge.op import nodal_min + return actx.to_numpy(nodal_min(dcoll, "vol", x))[()] + + def vol_max(x): + from grudge.op import nodal_max + return actx.to_numpy(nodal_max(dcoll, "vol", x))[()] + + def my_write_status(cv, dv, dt, cfl): + status_msg = f"-------- dt = {dt:1.3e}, cfl = {cfl:1.4f}" + p_min = vol_min(dv.pressure) + p_max = vol_max(dv.pressure) + t_min = vol_min(dv.temperature) + t_max = vol_max(dv.temperature) + + dv_status_msg = ( + f"\n-------- P (min, max) (Pa) = ({p_min:1.9e}, {p_max:1.9e})") + dv_status_msg += ( + f"\n-------- T (min, max) (K) = ({t_min:7g}, {t_max:7g})") + + status_msg += dv_status_msg + status_msg += "\n" + + if rank == 0: + logger.info(status_msg) + + from mirgecom.viscous import get_viscous_timestep, get_viscous_cfl + + def my_get_timestep(t, dt, state): + t_remaining = max(0, t_final - t) + if constant_cfl: + ts_field = current_cfl * get_viscous_timestep(dcoll, state=state) + from grudge.op import nodal_min + dt = actx.to_numpy(nodal_min(dcoll, "vol", ts_field)) + cfl = current_cfl + else: + ts_field = get_viscous_cfl(dcoll, dt=dt, state=state) + from grudge.op import nodal_max + cfl = actx.to_numpy(nodal_max(dcoll, "vol", ts_field)) + + return ts_field, cfl, min(t_remaining, dt) + + def my_write_viz(step, t, fluid_state): + cv = fluid_state.cv + dv = fluid_state.dv + mu = fluid_state.viscosity + + """ + exact_cv = initializer(x_vec=nodes, eos=gas_model.eos, time=t) + exact_smoothness = smoothness_indicator(dcoll, exact_cv.mass, + kappa=kappa, s0=s0) + exact_state = create_fluid_state(cv=exact_cv, + smoothness_mu=exact_smoothness) + + # try using the divergence to compute the smoothness field + #exact_grad_cv = grad_cv_operator_compiled(fluid_state=exact_state, + #time=t) + exact_grad_cv = grad_cv_operator(dcoll, gas_model, boundaries, + exact_state, time=current_t, + quadrature_tag=quadrature_tag) + from mirgecom.fluid import velocity_gradient + exact_grad_v = velocity_gradient(exact_cv, exact_grad_cv) + + # make a smoothness indicator + # try using the divergence to compute the smoothness field + exact_smoothness = compute_smoothness(exact_cv, exact_grad_cv) + + exact_state = create_fluid_state(cv=exact_cv, + smoothness_mu=exact_smoothness) + """ + + viz_fields = [("cv", cv), + ("dv", dv), + # ("exact_cv", exact_state.cv), + # ("exact_grad_v_x", exact_grad_v[0]), + # ("exact_grad_v_y", exact_grad_v[1]), + # ("exact_dv", exact_state.dv), + ("mu", mu)] + from mirgecom.simutil import write_visfile + write_visfile(dcoll, viz_fields, visualizer, vizname=vizname, + 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): + # Note: This health check is tuned s.t. it is a test that + # the case gets the expected solution. If dt,t_final or + # other run parameters are changed, this check should + # be changed accordingly. + 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, + health_pres_min, health_pres_max), + op="lor"): + health_error = True + 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, + health_temp_min, health_temp_max), + 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=})") + + return health_error + + def my_pre_step(step, t, dt, state): + 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) + do_status = check_step(step=step, interval=nstatus) + + if any([do_viz, do_restart, do_health, do_status, constant_cfl]): + fluid_state = create_fluid_state(cv=state, + smoothness_mu=no_smoothness) + if use_av > 1: + # recompute the dv to have the correct smoothness + # this is forcing a recompile, only do it at dump time + # not sure why the compiled version of grad_cv doesn't work + if do_viz: + # use the divergence to compute the smoothness field + force_evaluation(actx, t) + grad_cv = grad_cv_operator_compiled(fluid_state, time=t) + smoothness = compute_smoothness_compiled(state, grad_cv) + + # this works, but seems a lot of work, + # not sure if it's really faster + # avoids re-computing the temperature + from dataclasses import replace + new_dv = replace(fluid_state.dv, smoothness_mu=smoothness) + fluid_state = replace(fluid_state, dv=new_dv) + new_tv = gas_model.transport.transport_vars( + cv=state, dv=new_dv, eos=gas_model.eos) + fluid_state = replace(fluid_state, tv=new_tv) + + dv = fluid_state.dv + # if the time integrator didn't force_eval, do so now + if not force_eval: + fluid_state = force_evaluation(actx, fluid_state) + dv = fluid_state.dv + + ts_field, cfl, dt = my_get_timestep(t, dt, fluid_state) + + if do_health: + health_errors = \ + global_reduce(my_health_check(state, dv), op="lor") + + if health_errors: + if rank == 0: + logger.info("Fluid solution failed health check.") + raise MyRuntimeError("Failed simulation health check.") + + if do_status: + my_write_status(dt=dt, cfl=cfl, dv=dv, cv=state) + + if do_restart: + my_write_restart(step=step, t=t, state=state) + + if do_viz: + my_write_viz(step=step, t=t, fluid_state=fluid_state) + + if constant_cfl: + dt = get_sim_timestep(dcoll, fluid_state, t, dt, + current_cfl, t_final, constant_cfl) + + except MyRuntimeError: + if rank == 0: + logger.info("Errors detected; attempting graceful exit.") + my_write_viz(step=step, t=t, fluid_state=fluid_state) + 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) + 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, quadrature_tag=quadrature_tag) + ) + + def _my_rhs_av(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, quadrature_tag=quadrature_tag) + + av_laplacian_operator(dcoll, fluid_state=fluid_state, + boundaries=boundaries, + time=t, gas_model=gas_model, + alpha=alpha, s0=s0, kappa=kappa, + quadrature_tag=quadrature_tag) + ) + + def _my_rhs_phys_visc_av(t, state): + + smoothness = smoothness_indicator(dcoll, state.mass, + kappa=kappa, s0=s0) + fluid_state = make_fluid_state(cv=state, gas_model=gas_model, + smoothness_mu=smoothness) + + return ( + ns_operator(dcoll, state=fluid_state, time=t, + boundaries=boundaries, + gas_model=gas_model, quadrature_tag=quadrature_tag) + ) + + def _my_rhs_phys_visc_div_av(t, state): + + fluid_state = make_fluid_state(cv=state, gas_model=gas_model, + smoothness_mu=no_smoothness) + + # use the divergence to compute the smoothness field + grad_cv = grad_cv_operator(dcoll, gas_model, boundaries, fluid_state, + time=t, quadrature_tag=quadrature_tag) + smoothness = compute_smoothness(state, grad_cv) + + from dataclasses import replace + new_dv = replace(fluid_state.dv, smoothness_mu=smoothness) + fluid_state = replace(fluid_state, dv=new_dv) + new_tv = gas_model.transport.transport_vars( + cv=state, dv=new_dv, eos=gas_model.eos) + fluid_state = replace(fluid_state, tv=new_tv) + + return ( + ns_operator(dcoll, state=fluid_state, time=t, + boundaries=boundaries, + gas_model=gas_model, quadrature_tag=quadrature_tag, + grad_cv=grad_cv) + ) + + my_rhs = (_my_rhs if use_av == 0 else _my_rhs_av if use_av == 1 else + _my_rhs_phys_visc_div_av) + + 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) + + # Dump the final data + if rank == 0: + logger.info("Checkpointing final state ...") + + if use_av < 2: + current_state = create_fluid_state(cv=current_cv) + else: + current_state = create_fluid_state(cv=current_cv, + smoothness_mu=no_smoothness) + + # use the divergence to compute the smoothness field + current_grad_cv = grad_cv_operator_compiled(current_state, + time=current_t) + smoothness = compute_smoothness_compiled(current_cv, + current_grad_cv) + from dataclasses import replace + new_dv = replace(current_state.dv, smoothness_mu=smoothness) + current_state = replace(current_state, dv=new_dv) + + final_dv = current_state.dv + ts_field, cfl, dt = my_get_timestep(t=current_t, dt=current_dt, + state=current_state) + my_write_status(dt=dt, cfl=cfl, cv=current_cv, dv=final_dv) + my_write_viz(step=current_step, t=current_t, fluid_state=current_state) + 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 = "doublemach" + 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("--log", action="store_true", default=True, + help="turn on logging") + parser.add_argument("--leap", action="store_true", + help="use leap timestepper") + 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() + lazy = args.lazy + if args.profiling: + if lazy: + raise ValueError("Can't use lazy and profiling together.") + + from grudge.array_context import get_reasonable_array_context_class + actx_class = get_reasonable_array_context_class(lazy=lazy, + distributed=True) + + 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(use_logmgr=args.log, use_leap=args.leap, use_profiling=args.profiling, + use_overintegration=args.overintegration, lazy=lazy, + casename=casename, rst_filename=rst_filename, actx_class=actx_class) + +# vim: foldmethod=marker diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index 282db9a17..e197cadc0 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -832,8 +832,6 @@ def __init__(self, boundary_gradient_cv_func=None, # Returns the boundary value for grad(temperature) boundary_gradient_temperature_func=None, - # For artificial viscosity - grad fluid soln on boundary - boundary_grad_av_func=None, ): """Initialize the PrescribedFluidBoundary and methods.""" self._bnd_state_func = boundary_state_func @@ -1191,7 +1189,10 @@ def farfield_state(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): ) return make_fluid_state(cv=cv_infinity, gas_model=gas_model, - temperature_seed=free_stream_temperature) + temperature_seed=free_stream_temperature, + smoothness_mu=state_minus.smoothness_mu, + smoothness_kappa=state_minus.smoothness_kappa, + smoothness_beta=state_minus.smoothness_beta) def temperature_bc(self, state_minus, **kwargs): """Return farfield temperature for use in grad(temperature).""" @@ -1312,7 +1313,10 @@ def outflow_state(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): species_mass=state_minus.cv.species_mass) return make_fluid_state(cv=cv_outflow, gas_model=gas_model, - temperature_seed=state_minus.temperature) + temperature_seed=state_minus.temperature, + smoothness_mu=state_minus.smoothness_mu, + smoothness_kappa=state_minus.smoothness_kappa, + smoothness_beta=state_minus.smoothness_beta) def outflow_state_for_diffusion(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): @@ -1355,7 +1359,10 @@ def outflow_state_for_diffusion(self, dcoll, dd_bdry, gas_model, species_mass=state_minus.species_mass_density ) return make_fluid_state(cv=cv_plus, gas_model=gas_model, - temperature_seed=state_minus.temperature) + temperature_seed=state_minus.temperature, + smoothness_mu=state_minus.smoothness_mu, + smoothness_kappa=state_minus.smoothness_kappa, + smoothness_beta=state_minus.smoothness_beta) def inviscid_boundary_flux(self, dcoll, dd_bdry, gas_model, state_minus, numerical_flux_func=inviscid_facial_flux_rusanov, **kwargs): @@ -1497,7 +1504,10 @@ def inflow_state(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): species_mass=species_mass_boundary) return make_fluid_state(cv=boundary_cv, gas_model=gas_model, - temperature_seed=state_minus.temperature) + temperature_seed=state_minus.temperature, + smoothness_mu=state_minus.smoothness_mu, + smoothness_kappa=state_minus.smoothness_kappa, + smoothness_beta=state_minus.smoothness_beta) class RiemannOutflowBoundary(PrescribedFluidBoundary): @@ -1599,7 +1609,10 @@ def outflow_state(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): species_mass=species_mass_boundary) return make_fluid_state(cv=boundary_cv, gas_model=gas_model, - temperature_seed=state_minus.temperature) + temperature_seed=state_minus.temperature, + smoothness_mu=state_minus.smoothness_mu, + smoothness_kappa=state_minus.smoothness_kappa, + smoothness_beta=state_minus.smoothness_beta) class IsothermalWallBoundary(MengaldoBoundaryCondition): @@ -1825,4 +1838,7 @@ def outflow_state(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs): ) return make_fluid_state(cv=boundary_cv, gas_model=gas_model, - temperature_seed=state_minus.temperature) + temperature_seed=state_minus.temperature, + smoothness_mu=state_minus.smoothness_mu, + smoothness_kappa=state_minus.smoothness_kappa, + smoothness_beta=state_minus.smoothness_beta) diff --git a/mirgecom/eos.py b/mirgecom/eos.py index 6d8ba5625..032b89721 100644 --- a/mirgecom/eos.py +++ b/mirgecom/eos.py @@ -74,11 +74,18 @@ class GasDependentVars: .. attribute:: temperature .. attribute:: pressure + .. attribute:: speed_of_sound + .. attribute:: smoothness_mu + .. attribute:: smoothness_kappa + .. attribute:: smoothness_beta """ temperature: DOFArray pressure: DOFArray speed_of_sound: DOFArray + smoothness_mu: DOFArray + smoothness_kappa: DOFArray + smoothness_beta: DOFArray @dataclass_array_container @@ -154,7 +161,7 @@ def kinetic_energy(self, cv: ConservedVars): """Get the kinetic energy for the gas.""" @abstractmethod - def gamma(self, cv: ConservedVars, temperature: DOFArray): + def gamma(self, cv: ConservedVars, temperature: DOFArray = None): """Get the ratio of gas specific heats Cp/Cv.""" @abstractmethod @@ -167,7 +174,10 @@ def get_internal_energy(self, temperature, species_mass_fractions=None): def dependent_vars( self, cv: ConservedVars, - temperature_seed: Optional[DOFArray] = None) -> GasDependentVars: + temperature_seed: Optional[DOFArray] = None, + smoothness_mu: Optional[DOFArray] = None, + smoothness_kappa: Optional[DOFArray] = None, + smoothness_beta: Optional[DOFArray] = None) -> GasDependentVars: """Get an agglomerated array of the dependent variables. Certain implementations of :class:`GasEOS` (e.g. :class:`MixtureEOS`) @@ -175,10 +185,23 @@ def dependent_vars( given. """ temperature = self.temperature(cv, temperature_seed) + # MJA, it doesn't appear that we can have a None field embedded inside DV, + # make a dummy smoothness in this case + zeros = cv.array_context.zeros_like(cv.mass) + if smoothness_mu is None: + smoothness_mu = zeros + if smoothness_kappa is None: + smoothness_kappa = zeros + if smoothness_beta is None: + smoothness_beta = zeros + return GasDependentVars( temperature=temperature, pressure=self.pressure(cv, temperature), - speed_of_sound=self.sound_speed(cv, temperature) + speed_of_sound=self.sound_speed(cv, temperature), + smoothness_mu=smoothness_mu, + smoothness_kappa=smoothness_kappa, + smoothness_beta=smoothness_beta ) @@ -229,7 +252,10 @@ def get_species_source_terms(self, cv: ConservedVars, temperature: DOFArray): def dependent_vars( self, cv: ConservedVars, - temperature_seed: Optional[DOFArray] = None) -> MixtureDependentVars: + temperature_seed: Optional[DOFArray] = None, + smoothness_mu: Optional[DOFArray] = None, + smoothness_kappa: Optional[DOFArray] = None, + smoothness_beta: Optional[DOFArray] = None) -> MixtureDependentVars: """Get an agglomerated array of the dependent variables. Certain implementations of :class:`GasEOS` (e.g. :class:`MixtureEOS`) @@ -237,11 +263,24 @@ def dependent_vars( given. """ temperature = self.temperature(cv, temperature_seed) + # MJA, it doesn't appear that we can have a None field embedded inside DV, + # make a dummy smoothness in this case + zeros = cv.array_context.zeros_like(cv.mass) + if smoothness_mu is None: + smoothness_mu = zeros + if smoothness_kappa is None: + smoothness_kappa = zeros + if smoothness_beta is None: + smoothness_beta = zeros + return MixtureDependentVars( temperature=temperature, pressure=self.pressure(cv, temperature), speed_of_sound=self.sound_speed(cv, temperature), - species_enthalpies=self.species_enthalpies(cv, temperature) + species_enthalpies=self.species_enthalpies(cv, temperature), + smoothness_mu=smoothness_mu, + smoothness_kappa=smoothness_kappa, + smoothness_beta=smoothness_beta ) @@ -606,7 +645,7 @@ def heat_capacity_cv(self, cv: ConservedVars, temperature): / self.gamma(cv, temperature) ) - def gamma(self, cv: ConservedVars, temperature): + def gamma(self, cv: ConservedVars, temperature): # type: ignore[override] r"""Get mixture-averaged heat capacity ratio, $\frac{C_p}{C_p - R_s}$. Parameters @@ -726,7 +765,9 @@ def get_species_molecular_weights(self): def species_enthalpies(self, cv: ConservedVars, temperature): """Get the species specific enthalpies.""" - return self._pyrometheus_mech.get_species_enthalpies_rt(temperature) + spec_r = self._pyrometheus_mech.gas_constant/self._pyrometheus_mech.wts + return (spec_r * temperature + * self._pyrometheus_mech.get_species_enthalpies_rt(temperature)) def get_production_rates(self, cv: ConservedVars, temperature): r"""Get the production rate for each species. diff --git a/mirgecom/gas_model.py b/mirgecom/gas_model.py index 5782fe320..7ae828067 100644 --- a/mirgecom/gas_model.py +++ b/mirgecom/gas_model.py @@ -112,6 +112,9 @@ class FluidState: .. autoattribute:: nspecies .. autoattribute:: pressure .. autoattribute:: temperature + .. autoattribute:: smoothness_mu + .. autoattribute:: smoothness_kappa + .. autoattribute:: smoothness_beta .. autoattribute:: velocity .. autoattribute:: speed .. autoattribute:: wavespeed @@ -152,6 +155,21 @@ def temperature(self): """Return the gas temperature.""" return self.dv.temperature + @property + def smoothness_mu(self): + """Return the smoothness_mu field.""" + return self.dv.smoothness_mu + + @property + def smoothness_kappa(self): + """Return the smoothness_kappa field.""" + return self.dv.smoothness_kappa + + @property + def smoothness_beta(self): + """Return the smoothness_beta field.""" + return self.dv.smoothness_beta + @property def mass_density(self): """Return the gas density.""" @@ -258,8 +276,11 @@ def species_diffusivity(self): return self.tv.species_diffusivity -def make_fluid_state(cv, gas_model, temperature_seed=None, limiter_func=None, - limiter_dd=None): +def make_fluid_state(cv, gas_model, temperature_seed=None, + smoothness_mu=None, + smoothness_kappa=None, + smoothness_beta=None, + limiter_func=None, limiter_dd=None): """Create a fluid state from the conserved vars and physical gas model. Parameters @@ -290,15 +311,27 @@ def make_fluid_state(cv, gas_model, temperature_seed=None, limiter_func=None, """ temperature = gas_model.eos.temperature(cv, temperature_seed=temperature_seed) pressure = gas_model.eos.pressure(cv, temperature=temperature) + actx = cv.array_context if limiter_func: cv = limiter_func(cv=cv, pressure=pressure, temperature=temperature, dd=limiter_dd) + # FIXME work-around for now + smoothness_mu = (actx.zeros_like(pressure) if smoothness_mu + is None else smoothness_mu) + smoothness_kappa = (actx.zeros_like(pressure) if smoothness_kappa + is None else smoothness_kappa) + smoothness_beta = (actx.zeros_like(pressure) if smoothness_beta + is None else smoothness_beta) + dv = GasDependentVars( temperature=temperature, pressure=pressure, - speed_of_sound=gas_model.eos.sound_speed(cv, temperature) + speed_of_sound=gas_model.eos.sound_speed(cv, temperature), + smoothness_mu=smoothness_mu, + smoothness_kappa=smoothness_kappa, + smoothness_beta=smoothness_beta ) from mirgecom.eos import MixtureEOS, MixtureDependentVars @@ -307,11 +340,15 @@ def make_fluid_state(cv, gas_model, temperature_seed=None, limiter_func=None, temperature=dv.temperature, pressure=dv.pressure, speed_of_sound=dv.speed_of_sound, + smoothness_mu=dv.smoothness_mu, + smoothness_kappa=dv.smoothness_kappa, + smoothness_beta=dv.smoothness_beta, species_enthalpies=gas_model.eos.species_enthalpies(cv, temperature)) if gas_model.transport is not None: tv = gas_model.transport.transport_vars(cv=cv, dv=dv, eos=gas_model.eos) return ViscousFluidState(cv=cv, dv=dv, tv=tv) + return FluidState(cv=cv, dv=dv) @@ -363,8 +400,21 @@ 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) + smoothness_mu = None + if state.dv.smoothness_mu is not None: + smoothness_mu = op.project(dcoll, src, tgt, state.dv.smoothness_mu) + smoothness_kappa = None + if state.dv.smoothness_kappa is not None: + smoothness_kappa = op.project(dcoll, src, tgt, state.dv.smoothness_kappa) + smoothness_beta = None + if state.dv.smoothness_beta is not None: + smoothness_beta = op.project(dcoll, src, tgt, state.dv.smoothness_beta) + return make_fluid_state(cv=cv_sd, gas_model=gas_model, temperature_seed=temperature_seed, + smoothness_mu=smoothness_mu, + smoothness_kappa=smoothness_kappa, + smoothness_beta=smoothness_beta, limiter_func=limiter_func, limiter_dd=tgt) @@ -376,6 +426,9 @@ def _getattr_ish(obj, name): def make_fluid_state_trace_pairs(cv_pairs, gas_model, temperature_seed_pairs=None, + smoothness_mu_pairs=None, + smoothness_kappa_pairs=None, + smoothness_beta_pairs=None, limiter_func=None): """Create a fluid state from the conserved vars and equation of state. @@ -414,15 +467,37 @@ def make_fluid_state_trace_pairs(cv_pairs, gas_model, temperature_seed_pairs=Non from grudge.trace_pair import TracePair if temperature_seed_pairs is None: temperature_seed_pairs = [None] * len(cv_pairs) + if smoothness_mu_pairs is None: + smoothness_mu_pairs = [None] * len(cv_pairs) + if smoothness_kappa_pairs is None: + smoothness_kappa_pairs = [None] * len(cv_pairs) + if smoothness_beta_pairs is None: + smoothness_beta_pairs = [None] * len(cv_pairs) return [TracePair( cv_pair.dd, - interior=make_fluid_state(cv_pair.int, gas_model, - temperature_seed=_getattr_ish(tseed_pair, "int"), - limiter_func=limiter_func, limiter_dd=cv_pair.dd), - exterior=make_fluid_state(cv_pair.ext, gas_model, - temperature_seed=_getattr_ish(tseed_pair, "ext"), - limiter_func=limiter_func, limiter_dd=cv_pair.dd)) - for cv_pair, tseed_pair in zip(cv_pairs, temperature_seed_pairs)] + interior=make_fluid_state( + cv_pair.int, gas_model, + temperature_seed=_getattr_ish(tseed_pair, "int"), + smoothness_mu=_getattr_ish(smoothness_mu_pair, "int"), + smoothness_kappa=_getattr_ish(smoothness_kappa_pair, "int"), + smoothness_beta=_getattr_ish(smoothness_beta_pair, "int"), + limiter_func=limiter_func, limiter_dd=cv_pair.dd), + exterior=make_fluid_state( + cv_pair.ext, gas_model, + temperature_seed=_getattr_ish(tseed_pair, "ext"), + smoothness_mu=_getattr_ish(smoothness_mu_pair, "ext"), + smoothness_kappa=_getattr_ish(smoothness_kappa_pair, "ext"), + smoothness_beta=_getattr_ish(smoothness_beta_pair, "ext"), + limiter_func=limiter_func, limiter_dd=cv_pair.dd)) + for cv_pair, + tseed_pair, + smoothness_mu_pair, + smoothness_kappa_pair, + smoothness_beta_pair in zip(cv_pairs, + temperature_seed_pairs, + smoothness_mu_pairs, + smoothness_kappa_pairs, + smoothness_beta_pairs)] class _FluidCVTag: @@ -433,6 +508,18 @@ class _FluidTemperatureTag: pass +class _FluidSmoothnessMuTag: + pass + + +class _FluidSmoothnessKappaTag: + pass + + +class _FluidSmoothnessBetaTag: + pass + + 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): @@ -537,9 +624,34 @@ def make_operator_fluid_states( dcoll, volume_state.temperature, volume_dd=dd_vol, comm_tag=(_FluidTemperatureTag, comm_tag))] + smoothness_mu_interior_pairs = None + if volume_state.smoothness_mu is not None: + smoothness_mu_interior_pairs = [ + interp_to_surf_quad(tpair=tpair) + for tpair in interior_trace_pairs( + dcoll, volume_state.smoothness_mu, volume_dd=dd_vol, + tag=(_FluidSmoothnessMuTag, comm_tag))] + smoothness_kappa_interior_pairs = None + if volume_state.smoothness_kappa is not None: + smoothness_kappa_interior_pairs = [ + interp_to_surf_quad(tpair=tpair) + for tpair in interior_trace_pairs( + dcoll, volume_state.smoothness_kappa, volume_dd=dd_vol, + tag=(_FluidSmoothnessKappaTag, comm_tag))] + smoothness_beta_interior_pairs = None + if volume_state.smoothness_beta is not None: + smoothness_beta_interior_pairs = [ + interp_to_surf_quad(tpair=tpair) + for tpair in interior_trace_pairs( + dcoll, volume_state.smoothness_beta, volume_dd=dd_vol, + tag=(_FluidSmoothnessBetaTag, comm_tag))] + interior_boundary_states_quad = \ make_fluid_state_trace_pairs(cv_interior_pairs, gas_model, tseed_interior_pairs, + smoothness_mu_interior_pairs, + smoothness_kappa_interior_pairs, + smoothness_beta_interior_pairs, limiter_func=limiter_func) # Interpolate the fluid state to the volume quadrature grid @@ -626,5 +738,8 @@ def replace_fluid_state( cv=new_cv, gas_model=gas_model, temperature_seed=new_tseed, + smoothness_mu=state.smoothness_mu, + smoothness_kappa=state.smoothness_kappa, + smoothness_beta=state.smoothness_beta, limiter_func=limiter_func, limiter_dd=limiter_dd) diff --git a/mirgecom/initializers.py b/mirgecom/initializers.py index e6b235853..4513c5105 100644 --- a/mirgecom/initializers.py +++ b/mirgecom/initializers.py @@ -11,9 +11,10 @@ .. autoclass:: Uniform .. autoclass:: AcousticPulse .. autoclass:: MixtureInitializer -.. autoclass:: PlanarDiscontinuity .. autoclass:: PlanarPoiseuille .. autoclass:: ShearFlow +.. autoclass:: PlanarDiscontinuity +.. autoclass:: MulticomponentTrig State Initializers ^^^^^^^^^^^^^^^^^^ @@ -367,7 +368,7 @@ class DoubleMachReflection: """ def __init__( - self, shock_location=1.0/6.0, shock_speed=4.0 + self, shock_location=1.0/6.0, shock_speed=4.0, shock_sigma=0.05 ): """Initialize double shock reflection parameters. @@ -377,9 +378,12 @@ def __init__( initial location of shock shock_speed: float shock speed, Mach number + shock_sigma: float + initial condition sharpness """ self._shock_location = shock_location self._shock_speed = shock_speed + self._shock_sigma = shock_sigma def __call__(self, x_vec, *, eos=None, time=0, **kwargs): r""" @@ -431,7 +435,7 @@ def __call__(self, x_vec, *, eos=None, time=0, **kwargs): xinter = (self._shock_location + y_rel/np.sqrt(3.0) + 2.0*self._shock_speed*t/np.sqrt(3.0)) - sigma = 0.05 + sigma = self._shock_sigma xtanh = 1.0/sigma*(x_rel-xinter) mass = rhol/2.0*(actx.np.tanh(-xtanh)+1.0)+rhor/2.0*(actx.np.tanh(xtanh)+1.0) rhoe = (rhoel/2.0*(actx.np.tanh(-xtanh)+1.0) @@ -635,7 +639,7 @@ def __init__( rho0=1.0, p0=1.0, center=None, velocity=None, spec_y0s=None, spec_amplitudes=None, - spec_centers=None + spec_centers=None, sigma=1.0 ): r"""Initialize MulticomponentLump parameters. @@ -652,6 +656,8 @@ def __init__( velocity: numpy.ndarray fixed flow velocity used for exact solution at t != 0, shape ``(dim,)`` + sigma: float + std deviation of the gaussian """ if center is None: center = np.zeros(shape=(dim,)) @@ -686,6 +692,7 @@ def __init__( self._spec_y0s = spec_y0s self._spec_centers = spec_centers self._spec_amplitudes = spec_amplitudes + self._sigma = sigma def __call__(self, x_vec, *, eos=None, time=0, **kwargs): """ @@ -726,8 +733,8 @@ def __call__(self, x_vec, *, eos=None, time=0, **kwargs): for i in range(self._nspecies): lump_loc = self._spec_centers[i] + loc_update rel_pos = x_vec - lump_loc - r2 = np.dot(rel_pos, rel_pos) - expterm = self._spec_amplitudes[i] * actx.np.exp(-r2) + r2 = np.dot(rel_pos, rel_pos)/(self._sigma**2) + expterm = self._spec_amplitudes[i] * actx.np.exp(-0.5*r2) species_mass[i] = self._rho0 * (self._spec_y0s[i] + expterm) return make_conserved(dim=self._dim, mass=mass, energy=energy, @@ -766,14 +773,159 @@ def exact_rhs(self, dcoll, cv, time=0.0): for i in range(self._nspecies): lump_loc = self._spec_centers[i] + loc_update rel_pos = nodes - lump_loc - r2 = np.dot(rel_pos, rel_pos) - expterm = self._spec_amplitudes[i] * actx.np.exp(-r2) - specrhs[i] = 2 * self._rho0 * expterm * np.dot(rel_pos, v) + r2 = np.dot(rel_pos, rel_pos)/self._sigma**2 + expterm = self._spec_amplitudes[i] * actx.np.exp(-0.5*r2) + specrhs[i] = self._rho0 * expterm * np.dot(rel_pos, v) / self._sigma**2 return make_conserved(dim=self._dim, mass=massrhs, energy=energyrhs, momentum=momrhs, species_mass=specrhs) +class MulticomponentTrig: + r"""Initializer for trig-distributed species fractions. + + The trig lump is defined by: + + .. math:: + + \rho &= 1.0\\ + {\rho}\mathbf{V} &= {\rho}\mathbf{V}_0\\ + {\rho}E &= \frac{p_0}{(\gamma - 1)} + \frac{1}{2}\rho{|V_0|}^{2}\\ + {\rho~Y_\alpha} &= {\rho~Y_\alpha}_{0} + + {a_\alpha}\sin{\omega(\mathbf{r} - \mathbf{v}t)}, + + where $\mathbf{V}_0$ is the fixed velocity specified by the user at init time, + and $\gamma$ is taken from the equation-of-state object (eos). + + The user-specified vector of initial values (${{Y}_\alpha}_0$) + for the mass fraction of each species, *spec_y0s*, and $a_\alpha$ is the + user-specified vector of amplitudes for each species, *spec_amplitudes*, and + $c_\alpha$ is the user-specified origin for each species, *spec_centers*. + + A call to this object after creation/init creates the trig solution at a given + time (*t*) relative to the configured origin (*center*), wave_vector k, and + background flow velocity (*velocity*). + + .. automethod:: __init__ + .. automethod:: __call__ + """ + + def __init__( + self, *, dim=1, nspecies=0, + rho0=1.0, p0=1.0, gamma=1.4, + center=None, velocity=None, + spec_y0s=None, spec_amplitudes=None, + spec_centers=None, + spec_omegas=None, + spec_diffusivities=None, + wave_vector=None + ): + r"""Initialize MulticomponentLump parameters. + + Parameters + ---------- + dim: int + specify the number of dimensions for the lump + rho0: float + specifies the value of $\rho_0$ + p0: float + specifies the value of $p_0$ + center: numpy.ndarray + center of lump, shape ``(dim,)`` + velocity: numpy.ndarray + fixed flow velocity used for exact solution at t != 0, + shape ``(dim,)`` + """ + if center is None: + center = np.zeros(shape=(dim,)) + if velocity is None: + velocity = np.zeros(shape=(dim,)) + 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,)) + if spec_centers is None: + spec_centers = make_obj_array([np.zeros(shape=dim,) + for i in range(nspecies)]) + if spec_omegas is None: + spec_omegas = 2.*np.pi*np.ones(shape=(nspecies,)) + + if spec_amplitudes is None: + spec_amplitudes = np.ones(shape=(nspecies,)) + if spec_diffusivities is None: + spec_diffusivities = np.ones(shape=(nspecies,)) + + if wave_vector is None: + wave_vector = np.zeros(shape=(dim,)) + wave_vector[0] = 1 + + if len(spec_y0s) != nspecies or\ + len(spec_amplitudes) != nspecies or\ + len(spec_centers) != nspecies: + raise ValueError(f"Expected nspecies={nspecies} inputs.") + for i in range(nspecies): + if len(spec_centers[i]) != dim: + raise ValueError(f"Expected {dim}-dimensional " + f"inputs for spec_centers.") + + self._nspecies = nspecies + self._dim = dim + self._velocity = velocity + self._center = center + self._p0 = p0 + self._rho0 = rho0 + self._spec_y0s = spec_y0s + self._spec_centers = spec_centers + self._spec_amps = spec_amplitudes + self._gamma = gamma + self._spec_omegas = spec_omegas + self._d = spec_diffusivities + self._wave_vector = wave_vector + + def __call__(self, x_vec, *, eos=None, time=0, **kwargs): + """ + Create a multi-component lump solution at time *t* and locations *x_vec*. + + The solution at time *t* is created by advecting the component mass lump + at the user-specified constant, uniform velocity + (``MulticomponentLump._velocity``). + + Parameters + ---------- + time: float + Current time at which the solution is desired + x_vec: numpy.ndarray + Nodal coordinates + eos: :class:`mirgecom.eos.IdealSingleGas` + Equation of state class with method to supply gas *gamma*. + """ + t = time + if x_vec.shape != (self._dim,): + print(f"len(x_vec) = {len(x_vec)}") + print(f"self._dim = {self._dim}") + raise ValueError(f"Expected {self._dim}-dimensional inputs.") + # actx = x_vec[0].array_context + mass = 0 * x_vec[0] + self._rho0 + mom = self._velocity * mass + 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 + 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) + spec_y = self._spec_y0s[i] + self._spec_amps[i]*expterm*trigterm + spec_mass[i] = mass * spec_y + + return make_conserved(dim=self._dim, mass=mass, energy=energy, + momentum=mom, species_mass=spec_mass) + + class AcousticPulse: r"""Solution initializer for N-dimensional isentropic Gaussian acoustic pulse. @@ -1283,6 +1435,7 @@ def exact_grad(self, x_vec, eos, cv_exact): """Return the exact gradient of the Poiseuille state.""" y = x_vec[1] x = x_vec[0] + # FIXME: Symbolic infrastructure could perhaps do this better ones = x / x mass = cv_exact.mass diff --git a/mirgecom/transport.py b/mirgecom/transport.py index c2f982151..5d07e26ad 100644 --- a/mirgecom/transport.py +++ b/mirgecom/transport.py @@ -13,6 +13,8 @@ .. autoclass:: TransportModel .. autoclass:: SimpleTransport .. autoclass:: PowerLawTransport +.. autoclass:: ArtificialViscosityTransportDiv +.. autoclass:: ArtificialViscosityTransportDiv2 Exceptions ^^^^^^^^^^ @@ -95,17 +97,20 @@ class TransportModel: """ def bulk_viscosity(self, cv: ConservedVars, - dv: Optional[GasDependentVars] = None) -> DOFArray: + dv: Optional[GasDependentVars] = None, + eos: Optional[GasEOS] = None) -> DOFArray: r"""Get the bulk viscosity for the gas (${\mu}_{B}$).""" raise NotImplementedError() def viscosity(self, cv: ConservedVars, - dv: Optional[GasDependentVars] = None) -> DOFArray: + dv: Optional[GasDependentVars] = None, + eos: Optional[GasEOS] = None) -> DOFArray: r"""Get the gas dynamic viscosity, $\mu$.""" raise NotImplementedError() def volume_viscosity(self, cv: ConservedVars, - dv: Optional[GasDependentVars] = None) -> DOFArray: + dv: Optional[GasDependentVars] = None, + eos: Optional[GasEOS] = None) -> DOFArray: r"""Get the 2nd coefficent of viscosity, $\lambda$.""" raise NotImplementedError() @@ -126,8 +131,8 @@ def transport_vars(self, cv: ConservedVars, eos: Optional[GasEOS] = None) -> GasTransportVars: r"""Compute the transport properties from the conserved state.""" return GasTransportVars( - bulk_viscosity=self.bulk_viscosity(cv=cv, dv=dv), - viscosity=self.viscosity(cv=cv, dv=dv), + bulk_viscosity=self.bulk_viscosity(cv=cv, dv=dv, eos=eos), + viscosity=self.viscosity(cv=cv, dv=dv, eos=eos), thermal_conductivity=self.thermal_conductivity(cv=cv, dv=dv, eos=eos), species_diffusivity=self.species_diffusivity(cv=cv, dv=dv, eos=eos) ) @@ -157,17 +162,20 @@ def __init__(self, bulk_viscosity=0, viscosity=0, thermal_conductivity=0, self._d_alpha = species_diffusivity def bulk_viscosity(self, cv: ConservedVars, - dv: Optional[GasDependentVars] = None) -> DOFArray: + dv: Optional[GasDependentVars] = None, + eos: Optional[GasEOS] = None) -> DOFArray: r"""Get the bulk viscosity for the gas, $\mu_{B}$.""" return self._mu_bulk*(0*cv.mass + 1.0) def viscosity(self, cv: ConservedVars, - dv: Optional[GasDependentVars] = None) -> DOFArray: + dv: Optional[GasDependentVars] = None, + eos: Optional[GasEOS] = None) -> DOFArray: r"""Get the gas dynamic viscosity, $\mu$.""" return self._mu*(0*cv.mass + 1.0) def volume_viscosity(self, cv: ConservedVars, - dv: Optional[GasDependentVars] = None) -> DOFArray: + dv: Optional[GasDependentVars] = None, + eos: Optional[GasEOS] = None) -> DOFArray: r"""Get the 2nd viscosity coefficent, $\lambda$. In this transport model, the second coefficient of viscosity is defined as: @@ -238,7 +246,8 @@ def __init__(self, alpha=0.6, beta=4.093e-7, sigma=2.5, n=.666, self._lewis = lewis def bulk_viscosity(self, cv: ConservedVars, # type: ignore[override] - dv: GasDependentVars) -> DOFArray: + dv: GasDependentVars, + eos: Optional[GasEOS] = None) -> DOFArray: r"""Get the bulk viscosity for the gas, $\mu_{B}$. .. math:: @@ -250,7 +259,8 @@ def bulk_viscosity(self, cv: ConservedVars, # type: ignore[override] # TODO: Should this be memoized? Avoid multiple calls? def viscosity(self, cv: ConservedVars, # type: ignore[override] - dv: GasDependentVars) -> DOFArray: + dv: GasDependentVars, + eos: Optional[GasEOS] = None) -> DOFArray: r"""Get the gas dynamic viscosity, $\mu$. $\mu = \beta{T}^n$ @@ -258,7 +268,8 @@ def viscosity(self, cv: ConservedVars, # type: ignore[override] return self._beta * dv.temperature**self._n def volume_viscosity(self, cv: ConservedVars, # type: ignore[override] - dv: GasDependentVars) -> DOFArray: + dv: GasDependentVars, + eos: Optional[GasEOS] = None) -> DOFArray: r"""Get the 2nd viscosity coefficent, $\lambda$. In this transport model, the second coefficient of viscosity is defined as: @@ -277,7 +288,6 @@ def thermal_conductivity(self, cv: ConservedVars, # type: ignore[override] .. math:: \kappa = \sigma\mu{C}_{v} - """ return ( self._sigma * self.viscosity(cv, dv) @@ -304,3 +314,185 @@ def species_diffusivity(self, cv: ConservedVars, # type: ignore[override] cv.mass*self._lewis*eos.gamma(cv, dv.temperature)) ) return self._d_alpha*(0*cv.mass + 1.) + + +class ArtificialViscosityTransportDiv(TransportModel): + r"""Transport model for add artificial viscosity. + + Inherits from (and implements) :class:`TransportModel`. + + Takes a physical transport model and adds the artificial viscosity + contribution to it. Defaults to simple transport with inviscid settings. + This is equivalent to inviscid flow with artifical viscosity enabled. + + .. automethod:: __init__ + .. automethod:: bulk_viscosity + .. automethod:: viscosity + .. automethod:: volume_viscosity + .. automethod:: species_diffusivity + .. automethod:: thermal_conductivity + """ + + def __init__(self, + av_mu, av_prandtl, physical_transport=None, + av_species_diffusivity=None): + """Initialize uniform, constant transport properties.""" + if physical_transport is None: + self._physical_transport = SimpleTransport() + else: + self._physical_transport = physical_transport + + if av_species_diffusivity is None: + av_species_diffusivity = np.empty((0,), dtype=object) + + self._av_mu = av_mu + self._av_prandtl = av_prandtl + + def av_viscosity(self, cv, dv, eos): + r"""Get the artificial viscosity for the gas.""" + actx = cv.array_context + return self._av_mu*actx.np.sqrt(np.dot(cv.velocity, cv.velocity) + + dv.speed_of_sound**2) + + def bulk_viscosity(self, cv: ConservedVars, # type: ignore[override] + dv: GasDependentVars, + eos: GasEOS) -> DOFArray: + r"""Get the bulk viscosity for the gas, $\mu_{B}$.""" + return self._physical_transport.bulk_viscosity(cv, dv) + + def viscosity(self, cv: ConservedVars, # type: ignore[override] + dv: GasDependentVars, + eos: GasEOS) -> DOFArray: + r"""Get the gas dynamic viscosity, $\mu$.""" + return (dv.smoothness_mu*self.av_viscosity(cv, dv, eos) + + self._physical_transport.viscosity(cv, dv)) + + def volume_viscosity(self, cv: ConservedVars, # type: ignore[override] + dv: GasDependentVars, + eos: GasEOS) -> DOFArray: + r"""Get the 2nd viscosity coefficent, $\lambda$. + + In this transport model, the second coefficient of viscosity is defined as: + + $\lambda = \left(\mu_{B} - \frac{2\mu}{3}\right)$ + """ + return (dv.smoothness_mu*self.av_viscosity(cv, dv, eos) + + self._physical_transport.volume_viscosity(cv, dv)) + + def thermal_conductivity(self, cv: ConservedVars, # type: ignore[override] + dv: GasDependentVars, + eos: GasEOS) -> DOFArray: + r"""Get the gas thermal_conductivity, $\kappa$.""" + mu = self.av_viscosity(cv, dv, eos) + av_kappa = (dv.smoothness_mu*mu + * eos.heat_capacity_cp(cv, dv.temperature)/self._av_prandtl) + return av_kappa + self._physical_transport.thermal_conductivity( + cv, dv, eos) + + def species_diffusivity(self, cv: ConservedVars, + dv: Optional[GasDependentVars] = None, + eos: Optional[GasEOS] = None) -> DOFArray: + r"""Get the vector of species diffusivities, ${d}_{\alpha}$.""" + return self._physical_transport.species_diffusivity(cv, dv, eos) + + +class ArtificialViscosityTransportDiv2(TransportModel): + r"""Transport model for add artificial viscosity. + + Inherits from (and implements) :class:`TransportModel`. + + Takes a physical transport model and adds the artificial viscosity + contribution to it. Defaults to simple transport with inviscid settings. + This is equivalent to inviscid flow with artifical viscosity enabled. + + .. automethod:: __init__ + .. automethod:: bulk_viscosity + .. automethod:: viscosity + .. automethod:: volume_viscosity + .. automethod:: species_diffusivity + .. automethod:: thermal_conductivity + """ + + def __init__(self, + av_mu, av_kappa, av_beta, av_prandtl, + physical_transport=None, + av_species_diffusivity=None): + """Initialize uniform, constant transport properties.""" + if physical_transport is None: + self._physical_transport = SimpleTransport() + else: + self._physical_transport = physical_transport + + if av_species_diffusivity is None: + av_species_diffusivity = np.empty((0,), dtype=object) + + self._av_mu = av_mu + self._av_beta = av_beta + self._av_kappa = av_kappa + self._av_prandtl = av_prandtl + + def av_mu(self, cv, dv, eos): + r"""Get the shear artificial viscosity for the gas.""" + actx = cv.array_context + return (self._av_mu * cv.mass + * actx.np.sqrt(np.dot(cv.velocity, cv.velocity) + + dv.speed_of_sound**2)) + + def av_beta(self, cv, dv, eos): + r"""Get the shear artificial viscosity for the gas.""" + actx = cv.array_context + return (self._av_beta * cv.mass + * actx.np.sqrt(np.dot(cv.velocity, cv.velocity) + + dv.speed_of_sound**2)) + + def av_kappa(self, cv, dv, eos): + r"""Get the shear artificial viscosity for the gas.""" + actx = cv.array_context + return (self._av_kappa * cv.mass + * actx.np.sqrt(np.dot(cv.velocity, cv.velocity) + + dv.speed_of_sound**2)) + + def bulk_viscosity(self, cv: ConservedVars, # type: ignore[override] + dv: GasDependentVars, + eos: GasEOS) -> DOFArray: + r"""Get the bulk viscosity for the gas, $\mu_{B}$.""" + return (dv.smoothness_beta*self.av_beta(cv, dv, eos) + + self._physical_transport.bulk_viscosity(cv, dv)) + + def viscosity(self, cv: ConservedVars, # type: ignore[override] + dv: GasDependentVars, + eos: GasEOS) -> DOFArray: + r"""Get the gas dynamic viscosity, $\mu$.""" + return (dv.smoothness_mu*self.av_mu(cv, dv, eos) + + self._physical_transport.viscosity(cv, dv)) + + def volume_viscosity(self, cv: ConservedVars, # type: ignore[override] + dv: GasDependentVars, + eos: GasEOS) -> DOFArray: + r"""Get the 2nd viscosity coefficent, $\lambda$. + + In this transport model, the second coefficient of viscosity is: + + $\lambda = \left(\mu_{B} - \frac{2\mu}{3}\right)$ + """ + return (dv.smoothness_mu*self.av_mu(cv, dv, eos) + + self._physical_transport.volume_viscosity(cv, dv)) + + def thermal_conductivity(self, cv: ConservedVars, # type: ignore[override] + dv: GasDependentVars, + eos: GasEOS) -> DOFArray: + r"""Get the gas thermal_conductivity, $\kappa$.""" + cp = eos.heat_capacity_cp(cv, dv.temperature) + av_kappa = ( + cp*(dv.smoothness_beta*self.av_beta(cv, dv, eos)/self._av_prandtl + + dv.smoothness_kappa*self.av_kappa(cv, dv, eos)) + ) + + return (av_kappa + + self._physical_transport.thermal_conductivity(cv, dv, eos)) + + def species_diffusivity(self, cv: ConservedVars, + dv: Optional[GasDependentVars] = None, + eos: Optional[GasEOS] = None) -> DOFArray: + r"""Get the vector of species diffusivities, ${d}_{\alpha}$.""" + return self._physical_transport.species_diffusivity(cv, dv, eos) diff --git a/test/test_init.py b/test/test_init.py index c19dff5da..208a6b96d 100644 --- a/test/test_init.py +++ b/test/test_init.py @@ -391,7 +391,7 @@ def inf_norm(x): spec_r = make_obj_array([nodes - centers[i] for i in range(nspecies)]) r2 = make_obj_array([np.dot(spec_r[i], spec_r[i]) for i in range(nspecies)]) - expfactor = make_obj_array([spec_amplitudes[i] * actx.np.exp(- r2[i]) + expfactor = make_obj_array([spec_amplitudes[i] * actx.np.exp(-0.5*r2[i]) for i in range(nspecies)]) exp_mass = make_obj_array([rho0 * (spec_y0s[i] + expfactor[i]) for i in range(nspecies)])