Blast wave scaling test

Yurlungur commented 3 years ago

Problem setup

As descriged in the Athena test suite

https://www.astro.princeton.edu/~jstone/Athena/tests/blast/blast.html

and in

Zachary, Malagoli, A., & Colella,P., SIAM J. Sci. Comp., 15, 263 (1994); Balsara, D., & Spicer, D., JCP 149, 270 (1999); Londrillo, P. & Del Zanna, L., ApJ 530, 508 (2000).

Summary:

rho = 1
u = r > 0.1 ? 0.1*(3./2.) : 100.*0.1*(3./2.)
0 velocity
EOS = gamma law with gamma=5./3.

Scaling parameters

128^3 base grid on one node. Double resolution in one axis when doubling nodes.
4 levels of refinement
Meshblock sizes:
- 8^3
- 32^3
- 64^3
Problem setup:

Yurlungur commented 3 years ago

Refinement is 1st derivative of density, as set up by parthenon.

Not sure how long to run for. Suggestions welcome.

pgrete commented 3 years ago

Here's the version that matches Athena++ setup (used in the paper), i.e.,

128^3 root grid with 16^3 MeshBlocks
refinement on pressure gradient
run till t=0.1
note, that the ramp parameter is important to get the result reported in the paper. I'll ramp up the overpressure from r=0 to r=0.01

<comment>
problem   = spherical blast wave
reference = Gardiner. T.A. & Stone, J.M., JCP, 205, 509 (2005) (for MHD version of test)
configure = --prob=blast

<job>
problem_id = Blast      # problem ID: basename of output filenames

<output1>
file_type  = hst        # History data dump
dt         = 0.01       # time increment between outputs

<output2>
file_type  = hdf5       # Binary data dump
variable   = prim       # variables to be output
dt         = 0.10       # time increment between outputs

<time>
cfl_number = 0.3        # The Courant, Friedrichs, & Lewy (CFL) Number
nlim       = -1         # cycle limit
tlim       = 0.1        # time limit
integrator  = vl2       # time integration algorithm
xorder      = 2         # order of spatial reconstruction
ncycle_out  = 1         # interval for stdout summary info

<mesh>
nx1        = 128        # Number of zones in X1-direction
x1min      = -0.5       # minimum value of X1
x1max      = 0.5        # maximum value of X1
ix1_bc     = periodic   # inner-X1 boundary flag
ox1_bc     = periodic   # outer-X1 boundary flag

nx2        = 128        # Number of zones in X2-direction
x2min      = -0.5       # minimum value of X2
x2max      = 0.5        # maximum value of X2
ix2_bc     = periodic   # inner-X2 boundary flag
ox2_bc     = periodic   # outer-X2 boundary flag

nx3        = 128        # Number of zones in X3-direction
x3min      = -0.5       # minimum value of X3
x3max      = 0.5        # maximum value of X3
ix3_bc     = periodic   # inner-X3 boundary flag
ox3_bc     = periodic   # outer-X3 boundary flag

refinement  = adaptive # AMR
derefine_count = 10    # allow derefinement after 5 steps
numlevel    = 3        # number of AMR levels

<meshblock>
nx1        = 16        # Number of zones in X1-direction
nx2        = 16        # Number of zones in X2-direction
nx3        = 16        # Number of zones in X3-direction

<hydro>
gamma           = 1.666666666667 # gamma = C_p/C_v
iso_sound_speed = 0.4082482905   # equavalent to sqrt(gamma*p/d) for p=0.1, d=1

<problem>
compute_error = false  # check whether blast is spherical at end
pamb          = 0.001    # ambient pressure
prat          = 1.6e8  # Pressure ratio initially
radius        = 0.01    # Radius of the inner sphere
ramp          = 0.00
thr           = 0.1

And here's the corresponding input file for AthenaPK

# AthenaPK - a performance portable block structured AMR MHD code
# Copyright (c) 2020, Athena Parthenon Collaboration. All rights reserved.
# Licensed under the BSD 3-Clause License (the "LICENSE");

<comment>
problem = spherical blast wave

<job>
problem_id = blast

<parthenon/mesh>
refinement = adaptive
numlevel = 3

nx1 = 128
x1min = -0.5
x1max = 0.5
ix1_bc = periodic
ox1_bc = periodic

nx2 = 128
x2min = -0.5
x2max = 0.5
ix2_bc = periodic
ox2_bc = periodic

nx3 = 128
x3min = -0.5
x3max = 0.5
ix3_bc = periodic
ox3_bc = periodic

<parthenon/meshblock>
nx1=16
nx2=16
nx3=16

<parthenon/time>
integrator = vl2
cfl = 0.3
tlim = 0.10
nlim = 100000
perc_cycle_offset = 2 # number of inital cycles not to be included in perf calc
ncycle_out_mesh = -500

<hydro>
eos = adiabatic
riemann = hlle
reconstruction = plm
gamma = 1.666666666666667 # gamma = C_p/C_v
use_scratch = true 
scratch_level = 0 # 0 is actual scratch (tiny); 1 is HBM

<refinement>
type                         = pressure_gradient
threshold_pressure_gradient  = 0.1

<problem>
pamb       = 0.001      # ambient pressure
prat       = 1.6e8      # Pressure ratio initially
radius_outer  = 0.01    # Radius of the outer sphere
radius_inner  = 0.00    # Radius of the inter sphere (with ramp between inner and outer sphere)
damb       = 1.0        # ambient density
drat       = 1.0        # density ratio in sphere

<parthenon/output0>
file_type = hdf5
dt = 0.1 
variables = prim

pgrete commented 3 years ago

And here are some numbers: On 4x dual socket AMD EPYC 7H12 (i.e., 512 MPI ranks):

Athena++ 16^3 AMR: walltime 472 seconds at 2.38e8 zone-cycles/s
AthenaPK 16^3 AMR: walltime 692 seconds at 1.63e8 zone-cycles/s
AthenaPK 512^3 uni: walltime 3400 seconds at 4.22e8 zone-cycles/s

On 2x Power9 with 4x V100 each (i,.e., using 8x V100 total)

AthenaPK 512^3 uni: walltime 854 seconds at 16.8e8 zone-cycles/s
AthenaPK 16^3 AMR: tbd (but not great at the moment...)

Yurlungur commented 3 years ago

What is the pgen? Is it in the one in the Athena++ public version?

parthenon-hpc-lab / parthenon