|
| 1 | +__copyright__ = """ |
| 2 | +Copyright (C) 2021 University of Illinois Board of Trustees |
| 3 | +""" |
| 4 | + |
| 5 | +__license__ = """ |
| 6 | +Permission is hereby granted, free of charge, to any person obtaining a copy |
| 7 | +of this software and associated documentation files (the "Software"), to deal |
| 8 | +in the Software without restriction, including without limitation the rights |
| 9 | +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
| 10 | +copies of the Software, and to permit persons to whom the Software is |
| 11 | +furnished to do so, subject to the following conditions: |
| 12 | +
|
| 13 | +The above copyright notice and this permission notice shall be included in |
| 14 | +all copies or substantial portions of the Software. |
| 15 | +
|
| 16 | +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
| 17 | +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
| 18 | +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE |
| 19 | +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
| 20 | +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
| 21 | +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN |
| 22 | +THE SOFTWARE. |
| 23 | +""" |
| 24 | + |
| 25 | + |
| 26 | +from grudge.dof_desc import BoundaryDomainTag |
| 27 | +import pyopencl as cl |
| 28 | +import pyopencl.tools as cl_tools |
| 29 | + |
| 30 | +from arraycontext import thaw, freeze |
| 31 | +from grudge.array_context import PytatoPyOpenCLArrayContext |
| 32 | +from grudge.models.euler import ( |
| 33 | + EntropyStableEulerOperator, |
| 34 | + ConservedEulerField, |
| 35 | + PrescribedBC, |
| 36 | + conservative_to_primitive_vars, |
| 37 | +) |
| 38 | +from grudge.shortcuts import rk4_step |
| 39 | + |
| 40 | +from pytools.obj_array import make_obj_array |
| 41 | + |
| 42 | +import grudge.op as op |
| 43 | + |
| 44 | +import logging |
| 45 | +logger = logging.getLogger(__name__) |
| 46 | + |
| 47 | + |
| 48 | +def sod_shock_initial_condition(nodes, t=0): |
| 49 | + gamma = 1.4 |
| 50 | + dim = len(nodes) |
| 51 | + gmn1 = 1.0 / (gamma - 1.0) |
| 52 | + x = nodes[0] |
| 53 | + actx = x.array_context |
| 54 | + zeros = 0*x |
| 55 | + |
| 56 | + _x0 = 0.5 |
| 57 | + _rhoin = 1.0 |
| 58 | + _rhoout = 0.125 |
| 59 | + _pin = 1.0 |
| 60 | + _pout = 0.1 |
| 61 | + rhoin = zeros + _rhoin |
| 62 | + rhoout = zeros + _rhoout |
| 63 | + energyin = zeros + gmn1 * _pin |
| 64 | + energyout = zeros + gmn1 * _pout |
| 65 | + |
| 66 | + x0 = zeros + _x0 |
| 67 | + sigma = 1e-13 |
| 68 | + weight = 0.5 * (1.0 - actx.np.tanh(1.0/sigma * (x - x0))) |
| 69 | + |
| 70 | + mass = rhoout + (rhoin - rhoout)*weight |
| 71 | + energy = energyout + (energyin - energyout)*weight |
| 72 | + momentum = make_obj_array([zeros for _ in range(dim)]) |
| 73 | + |
| 74 | + return ConservedEulerField(mass=mass, energy=energy, momentum=momentum) |
| 75 | + |
| 76 | + |
| 77 | +def run_sod_shock_tube( |
| 78 | + actx, order=4, resolution=32, final_time=0.2, visualize=False): |
| 79 | + |
| 80 | + logger.info( |
| 81 | + """ |
| 82 | + Sod 1-D parameters:\n |
| 83 | + order: %s\n |
| 84 | + final_time: %s\n |
| 85 | + resolution: %s\n |
| 86 | + visualize: %s |
| 87 | + """, |
| 88 | + order, final_time, resolution, visualize |
| 89 | + ) |
| 90 | + |
| 91 | + # eos-related parameters |
| 92 | + gamma = 1.4 |
| 93 | + |
| 94 | + # {{{ discretization |
| 95 | + |
| 96 | + from meshmode.mesh.generation import generate_regular_rect_mesh |
| 97 | + |
| 98 | + dim = 1 |
| 99 | + box_ll = 0.0 |
| 100 | + box_ur = 1.0 |
| 101 | + mesh = generate_regular_rect_mesh( |
| 102 | + a=(box_ll,)*dim, |
| 103 | + b=(box_ur,)*dim, |
| 104 | + nelements_per_axis=(resolution,)*dim, |
| 105 | + boundary_tag_to_face={ |
| 106 | + "prescribed": ["+x", "-x"], |
| 107 | + } |
| 108 | + ) |
| 109 | + |
| 110 | + from grudge import DiscretizationCollection |
| 111 | + from grudge.dof_desc import \ |
| 112 | + DISCR_TAG_BASE, DISCR_TAG_QUAD |
| 113 | + from meshmode.discretization.poly_element import \ |
| 114 | + (default_simplex_group_factory, |
| 115 | + QuadratureSimplexGroupFactory) |
| 116 | + |
| 117 | + exp_name = f"fld-sod-1d-N{order}-K{resolution}" |
| 118 | + quad_tag = DISCR_TAG_QUAD |
| 119 | + |
| 120 | + dcoll = DiscretizationCollection( |
| 121 | + actx, mesh, |
| 122 | + discr_tag_to_group_factory={ |
| 123 | + DISCR_TAG_BASE: default_simplex_group_factory(dim, order), |
| 124 | + DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(order + 2) |
| 125 | + } |
| 126 | + ) |
| 127 | + |
| 128 | + # }}} |
| 129 | + |
| 130 | + # {{{ Euler operator |
| 131 | + |
| 132 | + dd_prescribed = BoundaryDomainTag("prescribed") |
| 133 | + bcs = { |
| 134 | + dd_prescribed: PrescribedBC(prescribed_state=sod_shock_initial_condition) |
| 135 | + } |
| 136 | + |
| 137 | + euler_operator = EntropyStableEulerOperator( |
| 138 | + dcoll, |
| 139 | + bdry_conditions=bcs, |
| 140 | + flux_type="lf", |
| 141 | + gamma=gamma, |
| 142 | + quadrature_tag=quad_tag |
| 143 | + ) |
| 144 | + |
| 145 | + def rhs(t, q): |
| 146 | + return euler_operator.operator(t, q) |
| 147 | + |
| 148 | + compiled_rhs = actx.compile(rhs) |
| 149 | + |
| 150 | + fields = sod_shock_initial_condition(thaw(dcoll.nodes(), actx)) |
| 151 | + |
| 152 | + from grudge.dt_utils import h_min_from_volume |
| 153 | + |
| 154 | + cfl = 0.01 |
| 155 | + cn = 0.5*(order + 1)**2 |
| 156 | + dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn |
| 157 | + |
| 158 | + logger.info("Timestep size: %g", dt) |
| 159 | + |
| 160 | + # }}} |
| 161 | + |
| 162 | + from grudge.shortcuts import make_visualizer |
| 163 | + |
| 164 | + vis = make_visualizer(dcoll) |
| 165 | + |
| 166 | + # {{{ time stepping |
| 167 | + |
| 168 | + step = 0 |
| 169 | + t = 0.0 |
| 170 | + while t < final_time: |
| 171 | + if step % 10 == 0: |
| 172 | + norm_q = actx.to_numpy(op.norm(dcoll, fields, 2)) |
| 173 | + logger.info("[%04d] t = %.5f |q| = %.5e", step, t, norm_q) |
| 174 | + if visualize: |
| 175 | + rho, velocity, pressure = \ |
| 176 | + conservative_to_primitive_vars(fields, gamma=gamma) |
| 177 | + vis.write_vtk_file( |
| 178 | + f"{exp_name}-{step:04d}.vtu", |
| 179 | + [ |
| 180 | + ("rho", rho), |
| 181 | + ("energy", fields.energy), |
| 182 | + ("momentum", fields.momentum), |
| 183 | + ("velocity", velocity), |
| 184 | + ("pressure", pressure) |
| 185 | + ] |
| 186 | + ) |
| 187 | + assert norm_q < 10000 |
| 188 | + |
| 189 | + fields = thaw(freeze(fields, actx), actx) |
| 190 | + fields = rk4_step(fields, t, dt, compiled_rhs) |
| 191 | + t += dt |
| 192 | + step += 1 |
| 193 | + |
| 194 | + # }}} |
| 195 | + |
| 196 | + |
| 197 | +def main(ctx_factory, order=4, final_time=0.2, resolution=32, visualize=False): |
| 198 | + cl_ctx = ctx_factory() |
| 199 | + queue = cl.CommandQueue(cl_ctx) |
| 200 | + actx = PytatoPyOpenCLArrayContext( |
| 201 | + queue, |
| 202 | + allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)), |
| 203 | + ) |
| 204 | + |
| 205 | + run_sod_shock_tube( |
| 206 | + actx, order=order, |
| 207 | + resolution=resolution, |
| 208 | + final_time=final_time, |
| 209 | + visualize=visualize) |
| 210 | + |
| 211 | + |
| 212 | +if __name__ == "__main__": |
| 213 | + import argparse |
| 214 | + |
| 215 | + parser = argparse.ArgumentParser() |
| 216 | + parser.add_argument("--order", default=4, type=int) |
| 217 | + parser.add_argument("--tfinal", default=0.2, type=float) |
| 218 | + parser.add_argument("--resolution", default=32, type=int) |
| 219 | + parser.add_argument("--visualize", action="store_true") |
| 220 | + args = parser.parse_args() |
| 221 | + |
| 222 | + logging.basicConfig(level=logging.INFO) |
| 223 | + main(cl.create_some_context, |
| 224 | + order=args.order, |
| 225 | + final_time=args.tfinal, |
| 226 | + resolution=args.resolution, |
| 227 | + visualize=args.visualize) |
0 commit comments