Basic Usage

python run.py path/to/config.yaml

When you run a simulation:

1. A timestamped folder is created in the same directory as the YAML file
2. The YAML config is copied into this run folder (original remains for easy re-running)
3. All outputs (GIFs, images, data) are saved to the run folder
4. Initial conditions (step 0) are always output along with subsequent steps

Example folder structure after a run:

unit_tests/
└── my_config_YYYY-MM-DD_HH:MM:SS/
    ├── my_config.yaml   # copied config
    └── <monitor outputs>   # GIFs, images, etc.

YAML Configuration Reference

Physics Mode

intsharp supports two physics modes:

• advection_only (default): Scalar advection with optional interface sharpening
• euler: Compressible Euler equations (1D single-phase/two-phase, plus 2D two-phase RTI in FV)

For Euler mode, use physics.mode: euler and configure material(s) and initial conditions. Single-phase example:

physics:
  mode: euler
  material:
    gamma: 1.4       # Ideal gas
    p_infinity: 0.0
  euler_initial_conditions:
    type: riemann
    x_discontinuity: 0.5
    left:  {rho: 1.0, u: 0.0, p: 1.0}
    right: {rho: 0.125, u: 0.0, p: 0.1}
  euler_bc: transmissive
  use_muscl: false   # or true for MUSCL + Barth-Jespersen

Two-phase (water-air) uses phase1/phase2 and two_phase_model: 5eq. 2D RTI uses directional BCs (euler_bc_x, euler_bc_y) and optional gravity (physics.gravity). See rayleigh_taylor_2d for the approved 2D example run.

physics:
  mode: euler
  phase1: {gamma: 1.4, p_infinity: 0.0, rho_ref: 2.0}
  phase2: {gamma: 1.4, p_infinity: 0.0, rho_ref: 1.0}
  two_phase_model: 5eq
  euler_spatial_discretization: fv
  euler_bc_x: periodic
  euler_bc_y: reflective
  gravity: {enabled: true, gx: 0.0, gy: -9.81}
  euler_initial_conditions:
    type: rti
    interface_y0: 0.5
    alpha1_top: 1.0
    alpha1_bottom: 0.0
    perturbation_amplitude: 0.01
    perturbation_mode_x: 1

Domain (1D)

1D domains are cell-centered: dx = L/n_points, with grid points at cell centers. This allows CFL = 1 to give exact return after one revolution (no phase drift).

domain:
  x_min: -0.5      # Left boundary
  x_max: 0.5       # Right boundary
  n_points: 100    # Number of grid points (dx = L/n = 0.01)

Domain (2D)

2D domains are cell-centered: dx = Lx/n_points_x, dy = Ly/n_points_y, with grid points at cell centers. This allows CFL = 1 in each direction to give exact return (no phase drift).

domain:
  x_min: -0.5        # Left boundary
  x_max: 0.5         # Right boundary
  n_points_x: 100    # Grid points in x (dx = Lx/n_points_x)
  y_min: -0.5        # Bottom boundary
  y_max: 0.5         # Top boundary
  n_points_y: 100    # Grid points in y (dy = Ly/n_points_y)

Time

time:
  dt: 0.01         # Time step size
  n_steps: 200     # Number of steps

Velocity

velocity: [0.5]        # 1D: advection velocity (list with 1 element)
velocity: [0.5, 0.0]   # 2D: (u, v) advection velocity

Fields

fields:
  - name: alpha
    initial_condition: "0.5 * (tanh((0.1 + x) / 0.02) + tanh((0.1 - x) / 0.02))"
    boundary:
      type: periodic   # or neumann, dirichlet
    sharpening: true           # optional per-field override
    sharpening_method: pm      # optional per-field method override

Per-field sharpening: Each field can override the global sharpening on/off via sharpening: true/false, and optionally select a different sharpening method via sharpening_method. This allows comparing multiple sharpening algorithms in the same simulation.

Expression-based Initial Conditions

Initial conditions can be Python expressions evaluated with:

• 1D: x (grid coordinates)
• 2D: x, y (meshgrid), r = sqrt(x² + y²) (radial distance)
• Math functions: sin, cos, tanh, exp, sqrt, pi, np (numpy)

Image-based Initial Conditions

Instead of an expression, you can specify an image file:

fields:
  - name: alpha
    initial_condition_image: "shapes/circle.png"
    boundary:
      type: periodic

The image is processed as follows:

1. Converted to grayscale
2. Resized to match the domain grid using nearest-neighbor interpolation
3. Thresholded: dark pixels (< 0.5) become α = 1, light pixels become α = 0
4. A tanh-smoothed interface is applied with ε = 3·dx

Image paths are resolved relative to the YAML config file location.

Boundary Condition Types

Type	Parameters
periodic	(none)
neumann	gradient_left, gradient_right
dirichlet	value_left, value_right

Solver

solver:
  type: upwind    # First-order upwind advection

For compressible runs, use:

physics:
  mode: euler
  euler_spatial_discretization: fv   # fv or dg
  flux_calculator: ausm_plus_up      # ausm_plus_up or hllc (fv)
  dg_order: 1                        # if dg: 1, 2, or 3
solver:
  type: upwind   # ignored in euler mode (kept for backward compatibility)

Timestepper

timestepper:
  type: euler     # or rk4

For Euler mode with euler_spatial_discretization: fv, this controls time integration directly.

For Euler mode with euler_spatial_discretization: dg, the current solver uses an internal SSP-RK3 update in intsharp/solvers/euler_dg_1d.py (the YAML timestepper is not used for DG at present).

Sharpening (Optional)

sharpening:
  enabled: true
  method: cl       # cl, pm, olsson_kreiss, acls, cls_2010,
                   # lcls_2012, lcls_2014, cls_2015, cls_2017, scls
  eps_target: 0.01 # Target interface thickness
  strength: 1.0    # Sharpening strength (Gamma)
  method_params:   # Optional method-specific parameters
    mapping_alpha: 2   # for cls_2010
    mapping_gamma: 2   # for cls_2015
    scls_alpha: 0.001  # for scls
    scls_beta: 1000    # for scls

Per-field method override: Each field can specify its own sharpening method via sharpening_method. This enables comparing multiple methods in a single simulation:

fields:
  - name: alpha_cl
    sharpening: true
    sharpening_method: cl
  - name: alpha_pm
    sharpening: true
    sharpening_method: pm

Output Monitors

When you use run.py, output is written to the timestamped run folder (the output.directory in the YAML is overridden). You can omit output.directory from your YAML.

output:
  monitors:
    - type: console      # Progress bar (no file output)
    - type: gif
      every_n_steps: 20
      field: alpha
    - type: gif          # Multi-field comparison
      every_n_steps: 20
      compare_fields:
        - field: alpha_no_sharp
          color: blue
          linestyle: "-"
        - field: alpha_pm
          color: red
          linestyle: "--"
    - type: mp4
      every_n_steps: 20
      field: alpha
    - type: png
      every_n_steps: 50
      field: alpha
    - type: svg
      every_n_steps: 50
      field: alpha
    - type: hdf5
      every_n_steps: 10
      fields: [alpha]

Monitor Types

Type	Output	Parameters
console	Progress bar (tqdm)	(none)
gif	Animated GIF	field or compare_fields, style, contour_levels
mp4	MP4 video	Same as gif (requires ffmpeg)
png	PNG images	field, every_n_steps or at_times
pdf	PDF images	field, every_n_steps or at_times
svg	SVG images	field, every_n_steps or at_times
hdf5	HDF5 file	fields (list), every_n_steps or at_times
txt	Columnar text	field, every_n_steps or at_times
curve	Curve format (x, y)	field, every_n_steps or at_times
metrics	metrics/*.tsv per field (step, t, ε_char, ε_char/ε_target, δ*_2, δ*_∞)	fields, optional interface_radius, metrics_eps_target, advection_velocity, initial_center

Advection runs with per-field sharpening convergence (convergence_tol set) also write convergence_times.tsv in the run directory: one row per sharpening field with t_f = simulation time when that field first satisfies the stop rule.

For 1D advection runs, the runner also exports interface_profiles/<field>_initial.tsv and interface_profiles/<field>_final.tsv (columns x, alpha) at step 0 and at the end of the run, for plotting initial vs. final interface shapes.

Advection-aware metrics: When the hat profile crosses the periodic boundary during advection, set advection_velocity on the metrics monitor so that the field is np.roll-ed back to center before computing ε_char and shape errors:

- type: metrics
  every_n_steps: 100
  fields: [alpha_cl]
  interface_radius: 0.025
  advection_velocity: 1.0     # constant velocity for periodic shift-back
  initial_center: 0.0         # optional, default 0.0

See gif for details on single-field vs multi-field modes and 2D visualization options.

Parameter Sweep & Stability Analysis

IntSharp includes tools for systematically evaluating sharpening methods across parameter space.

Running a sweep

# Full sweep (all 11 methods, 25x25 grid, n_substeps=[1,5,10])
make sweep

# Custom sweep
python scripts/sharpening_sweep.py --workers 8 --n-eps 30 --n-str 30 \
    --methods pm cl olsson_kreiss --substeps 1 10

Results are saved as CSV files in results/sweep/<method>/sweep_nsub<N>.csv.

Generating plots

make plot

This produces three heatmaps per method per n_substeps value:

• Stability map — green/red binary showing where the method is stable and bounded
• Sharpness map — transition cell count (lower = sharper), masked where unstable
• Conservation map — log₁₀(mass drift), masked where unstable

A summary grid PNG comparing all methods side-by-side is also generated.

Proposing new methods (LLM pipeline)

# Requires OPENAI_API_KEY environment variable
make propose

See sharpening_grammar.md for the allowed building blocks.