|
| 1 | +import jax.numpy as jnp |
| 2 | + |
| 3 | +# TODO: REMOVE THE FOLLOWING LINES |
| 4 | +import sys |
| 5 | + |
| 6 | +sys.path.append("../../") |
| 7 | + |
| 8 | +import adirondax as adx |
| 9 | +import time |
| 10 | +import matplotlib.pyplot as plt |
| 11 | + |
| 12 | +""" |
| 13 | +Simulate the Gresho Vortex |
| 14 | +
|
| 15 | +Philip Mocz (2025) |
| 16 | +""" |
| 17 | + |
| 18 | + |
| 19 | +def set_up_simulation(): |
| 20 | + # Define the parameters for the simulation |
| 21 | + nx = 128 |
| 22 | + nt = -1 |
| 23 | + t_stop = 3.0 |
| 24 | + |
| 25 | + params = { |
| 26 | + "physics": { |
| 27 | + "hydro": True, |
| 28 | + }, |
| 29 | + "mesh": { |
| 30 | + "type": "cartesian", |
| 31 | + "resolution": [nx, nx], |
| 32 | + "box_size": [1.0, 1.0], |
| 33 | + }, |
| 34 | + "time": { |
| 35 | + "span": t_stop, |
| 36 | + "num_timesteps": nt, |
| 37 | + }, |
| 38 | + "output": { |
| 39 | + "num_checkpoints": 100, |
| 40 | + "save": False, |
| 41 | + "plot_dynamic_range": 2.0, |
| 42 | + }, |
| 43 | + "hydro": { |
| 44 | + "eos": {"type": "ideal", "gamma": 5.0 / 3.0}, |
| 45 | + "slope_limiting": True, |
| 46 | + }, |
| 47 | + } |
| 48 | + |
| 49 | + # Initialize the simulation |
| 50 | + sim = adx.Simulation(params) |
| 51 | + |
| 52 | + # Set initial conditions |
| 53 | + sim.state["t"] = 0.0 |
| 54 | + X, Y = sim.mesh |
| 55 | + R = jnp.sqrt((X - 0.5) ** 2 + (Y - 0.5) ** 2) |
| 56 | + def v_phi(r): |
| 57 | + return jnp.where( |
| 58 | + r < 0.2, |
| 59 | + 5.0 * r, |
| 60 | + jnp.where(r < 0.4, 2.0 - 5.0 * r, 0.0), |
| 61 | + ) |
| 62 | + def P0(r): |
| 63 | + return jnp.where( |
| 64 | + r < 0.2, |
| 65 | + 5.0 + 12.5 * r**2, |
| 66 | + jnp.where(r < 0.4, 9.0 + 12.5 * r**2 - 20.0 * r + 4.0 * jnp.log(r / 0.2), 3.0 + 4.0 * jnp.log(2.0)), |
| 67 | + ) |
| 68 | + vx = -v_phi(R) * (Y - 0.5) / (R + (R==0)) |
| 69 | + vy = v_phi(R) * (X - 0.5) / (R + (R==0)) |
| 70 | + sim.state["rho"] = jnp.ones(X.shape) |
| 71 | + sim.state["vx"] = vx |
| 72 | + sim.state["vy"] = vy |
| 73 | + sim.state["P"] = P0(R) |
| 74 | + |
| 75 | + return sim |
| 76 | + |
| 77 | + |
| 78 | +def make_plot(sim): |
| 79 | + # Plot the solution |
| 80 | + plt.figure(figsize=(6, 4), dpi=80) |
| 81 | + v_phi = jnp.sqrt(sim.state["vx"]**2 + sim.state["vy"]**2) |
| 82 | + plt.imshow(v_phi.T, cmap="jet", vmin=0.0, vmax=1.0) |
| 83 | + plt.gca().invert_yaxis() |
| 84 | + plt.colorbar(label="v_phi") |
| 85 | + plt.tight_layout() |
| 86 | + plt.savefig("output.png", dpi=240) |
| 87 | + plt.show() |
| 88 | + |
| 89 | + |
| 90 | +def main(): |
| 91 | + sim = set_up_simulation() |
| 92 | + |
| 93 | + # Evolve the system |
| 94 | + t0 = time.time() |
| 95 | + sim.run() |
| 96 | + print("Run time (s): ", time.time() - t0) |
| 97 | + print("Steps taken:", sim.steps_taken) |
| 98 | + |
| 99 | + make_plot(sim) |
| 100 | + |
| 101 | + |
| 102 | +if __name__ == "__main__": |
| 103 | + main() |
0 commit comments