HMC Utilities - lattice/quda GitHub Wiki

#!/usr/bin/env bash
# run_dh_scaling.sh — dH-vs-dt scaling study on 8^4: 4 actions x 4 integrators.
# For each cell, runs HMC.Production once per step count (1 trajectory, no
# thermalization, fixed momentum seed => identical (U0, p0, phi) across the
# dt scan), extracts dH, and fits p in |dH| ~ dt^p by least squares on logs.
# Run from build/tests:  bash run_dh_scaling.sh
set -u

BIN=./hmc_test
RES_DIR="$PWD/tunecache_dh_scaling"
LOG_DIR="$PWD/logs_scaling"
DIM="--dim 8 8 8 8"
TAU=0.5
SEED=12345
TOL=1e-9
CELL_TIMEOUT=3600

[ -x "$BIN" ] || { echo "ERROR: run from build/tests"; exit 1; }
mkdir -p "$RES_DIR" "$LOG_DIR"
if [ -f "$RES_DIR/tunecache.tsv" ] && [ "$BIN" -nt "$RES_DIR/tunecache.tsv" ]; then
  echo "Wiping stale tunecache."
  rm -f "$RES_DIR"/*.tsv
fi

PREC_DBL="--prec double --prec-sloppy double --prec-precondition double"
PREC_MIX="--prec double --prec-sloppy single --prec-precondition single"   # nested (all actions)

action_flags() { # $1 action
  case "$1" in
    wilson)         echo "--dslash-type wilson --mass 0.5" ;;
    clover)         echo "--dslash-type clover --compute-clover true --clover-csw 1.0 --mass 0.5" ;;
    twisted-mass)   echo "--dslash-type twisted-mass --mu 0.1 --flavor singlet --mass 0.5" ;;
    twisted-clover) echo "--dslash-type twisted-clover --mu 0.1 --flavor singlet --compute-clover true --clover-csw 1.0 --mass 0.5" ;;
  esac
}
integ_id()    { case "$1" in leapfrog) echo 0;; omelyan) echo 1;; fgi) echo 2;; nested-fgi) echo 3;; esac; }
integ_steps() { case "$1" in leapfrog|omelyan) echo "16 24 36 54";; fgi|nested-fgi) echo "4 6 9 14";; esac; }
integ_p()     { case "$1" in leapfrog|omelyan) echo 2.0;; fgi|nested-fgi) echo 4.0;; esac; }

RESULTS=""

for action in wilson clover twisted-mass twisted-clover; do
  for integ in leapfrog omelyan fgi nested-fgi; do
    cell="${action}_${integ}"
    id=$(integ_id "$integ")
    expected=$(integ_p "$integ")
    note=""

    if [ "$integ" = nested-fgi ]; then
      mg="--mg-levels 2 --mg-block-size 0 4 4 4 4"
      prec=$PREC_MIX
    else
      mg="--mg-levels 1"
      prec=$PREC_DBL
    fi

    echo "=== [$(date '+%H:%M:%S')] $cell (expect p ~ $expected) $note"
    datafile="$LOG_DIR/$cell.dat"
    : > "$datafile"
    fail=0
    for n in $(integ_steps "$integ"); do
      log="$LOG_DIR/${cell}_n${n}.log"
      # shellcheck disable=SC2086
      QUDA_RESOURCE_PATH="$RES_DIR" timeout "$CELL_TIMEOUT" mpiexec -n 1 "$BIN" $DIM \
        $(action_flags "$action") $prec --tol $TOL \
        --hmc-integrator "$id" --hmc-n-steps "$n" --hmc-tau $TAU \
        --hmc-n-trajectories 1 --hmc-thermalization 0 --hmc-momentum-seed $SEED \
        $mg --gtest_filter=HMC.Production > "$log" 2>&1
      rc=$?
      dh=$(grep -oE 'dH = [+-][0-9.e+-]+' "$log" | head -1 | awk '{print $3}')
      if [ $rc -ne 0 ] || [ -z "$dh" ]; then
        echo "    n=$n : FAILED (rc=$rc)"
        fail=1
        continue
      fi
      dt=$(awk -v t=$TAU -v n=$n 'BEGIN{printf "%.8f", t/n}')
      echo "$dt $dh" >> "$datafile"
      printf '    n=%-3d dt=%-9s dH=%s\n' "$n" "$dt" "$dh"
    done

    npts=$(wc -l < "$datafile")
    if [ "$npts" -lt 3 ]; then
      RESULTS="$RESULTS$cell|FAIL|$expected|n/a|$note
"
      continue
    fi
    # least-squares fit of log|dH| vs log dt
    pfit=$(awk '{x=log($1); v=$2; if (v<0) v=-v; if (v<=0) next; y=log(v);
                 n++; sx+=x; sy+=y; sxx+=x*x; sxy+=x*y}
                END {if (n>=3) printf "%.3f", (n*sxy-sx*sy)/(n*sxx-sx*sx); else printf "n/a"}' "$datafile")
    verdict="PASS"
    awk -v p="$pfit" -v e="$expected" 'BEGIN{d=p-e; if (d<0) d=-d; exit !(d>0.5)}' && verdict="CHECK"
    [ $fail -eq 1 ] && verdict="PARTIAL"
    printf '    p_fit = %s (expected %s) -> %s\n' "$pfit" "$expected" "$verdict"
    RESULTS="$RESULTS$cell|$verdict|$expected|$pfit|$note
"
  done
done

{
  echo
  echo "==================== dH SCALING SUMMARY (8^4, tau=$TAU, tol=$TOL) ===================="
  printf '%-30s %-9s %-10s %-8s %s\n' "cell" "verdict" "expected p" "p_fit" "note"
  printf '%s\n' "--------------------------------------------------------------------------------"
  printf '%s' "$RESULTS" | while IFS='|' read -r c v e p nt; do
    [ -n "$c" ] && printf '%-30s %-9s %-10s %-8s %s\n' "$c" "$v" "$e" "$p" "$nt"
  done
  echo
  echo "CHECK = |p_fit - expected| > 0.5; inspect logs_scaling/<cell>.dat and per-n logs."
  echo "nested-fgi: effective 4th order on this step grid (h^4-dominated side of its"
  echo "two-term error crossover dH = a*(h/m)^2 + b*h^4; see run_nested_fgi_theory.sh)."
} | tee "$LOG_DIR/summary.txt"

#!/usr/bin/env bash
# run_dh_scaling_thermalized.sh — equilibrated dH-vs-dt scaling on 8^4.
#
# Protocol (per action x integrator x step size):
#   1. Load the 100-sweep quenched heatbath lattice (beta = 6.0).
#   2. 5 trajectories with forced accept   (dynamical-fermion thermalization)
#   3. 5 trajectories with Metropolis on   (equilibration, discarded)
#   4. 5 trajectories measured             (<|dH|> enters the scaling fit)
# The gauge is re-loaded fresh for every step size, so each dt point runs the
# identical protocol from the identical partially-thermalized, unbiased start.
# Fixed momentum seed => same heatbath momenta per trajectory index across the
# dt scan (variance reduction for the fit).
#
# Run from build/tests:  bash run_dh_scaling_thermalized.sh
set -u

BIN=./hmc_test
GAUGE="$PWD/therm_8x8_b6.0_hb100.lime"
RES_DIR="$PWD/tunecache_dh_scaling"
LOG_DIR="$PWD/logs_scaling_therm"
DIM="--dim 8 8 8 8"
TAU=0.5
SEED=12345
TOL=1e-9
N_TRAJ=15
N_THERM=5      # forced-accept trajectories
N_EQUIL=5      # Metropolis-on, discarded
N_SKIP=$((N_THERM + N_EQUIL))
CELL_TIMEOUT=3600

[ -x "$BIN" ] || { echo "ERROR: run from build/tests"; exit 1; }
[ -f "$GAUGE" ] || { echo "ERROR: thermalized gauge $GAUGE not found"; exit 1; }
mkdir -p "$RES_DIR" "$LOG_DIR"
if [ -f "$RES_DIR/tunecache.tsv" ] && [ "$BIN" -nt "$RES_DIR/tunecache.tsv" ]; then
  echo "Wiping stale tunecache."
  rm -f "$RES_DIR"/*.tsv
fi

PREC_DBL="--prec double --prec-sloppy double --prec-precondition double"
PREC_MIX="--prec double --prec-sloppy single --prec-precondition single"
PREC_SGL="--prec single --prec-sloppy single --prec-precondition single"

action_flags() {
  case "$1" in
    wilson)         echo "--dslash-type wilson --mass 0.5" ;;
    clover)         echo "--dslash-type clover --compute-clover true --clover-csw 1.0 --mass 0.5" ;;
    twisted-mass)   echo "--dslash-type twisted-mass --mu 0.1 --flavor singlet --mass 0.5" ;;
    twisted-clover) echo "--dslash-type twisted-clover --mu 0.1 --flavor singlet --compute-clover true --clover-csw 1.0 --mass 0.5" ;;
  esac
}
integ_id()    { case "$1" in leapfrog) echo 0;; omelyan) echo 1;; fgi) echo 2;; nested-fgi) echo 3;; esac; }
integ_steps() { case "$1" in leapfrog|omelyan) echo "16 24 36 54";; fgi|nested-fgi) echo "4 6 9 14";; esac; }
integ_p()     { case "$1" in leapfrog|omelyan) echo 2.0;; fgi|nested-fgi) echo 4.0;; esac; }

RESULTS=""

for action in wilson clover twisted-mass twisted-clover; do
  for integ in leapfrog omelyan fgi nested-fgi; do
    cell="${action}_${integ}"
    id=$(integ_id "$integ")
    expected=$(integ_p "$integ")
    note=""

    if [ "$integ" = nested-fgi ]; then
      mg="--mg-levels 2 --mg-block-size 0 4 4 4 4"
      # d/s/s mixed precision (in-place sloppy-clover family refresh keeps
      # clover-type actions at double-quality conservation)
      prec=$PREC_MIX
    else
      mg="--mg-levels 1"
      prec=$PREC_DBL
    fi

    echo "=== [$(date '+%H:%M:%S')] $cell (expect p ~ $expected) $note"
    datafile="$LOG_DIR/$cell.dat"
    : > "$datafile"
    fail=0
    for n in $(integ_steps "$integ"); do
      log="$LOG_DIR/${cell}_n${n}.log"
      # shellcheck disable=SC2086
      QUDA_RESOURCE_PATH="$RES_DIR" timeout "$CELL_TIMEOUT" mpiexec -n 1 "$BIN" $DIM \
        --load-gauge "$GAUGE" \
        $(action_flags "$action") $prec --tol $TOL \
        --hmc-integrator "$id" --hmc-n-steps "$n" --hmc-tau $TAU \
        --hmc-n-trajectories $N_TRAJ --hmc-thermalization $N_THERM --hmc-momentum-seed $SEED \
        $mg --gtest_filter=HMC.Production > "$log" 2>&1
      rc=$?
      # mean |dH| over the final measurement trajectories (skip therm + equil)
      stats=$(grep -E 'hmcRunQuda: (ACCEPTED|REJECTED)' "$log" \
        | awk -v skip=$N_SKIP 'NR>skip {v=$5+0; if (v<0) v=-v; s+=v; n++; if ($2=="ACCEPTED") a++}
                               END {if (n) printf "%.6e %d/%d", s/n, a, n; else printf "n/a n/a"}')
      dh=$(echo "$stats" | awk '{print $1}')
      acc=$(echo "$stats" | awk '{print $2}')
      if [ $rc -ne 0 ] || [ "$dh" = "n/a" ]; then
        echo "    n=$n : FAILED (rc=$rc)"
        fail=1
        continue
      fi
      dt=$(awk -v t=$TAU -v n=$n 'BEGIN{printf "%.8f", t/n}')
      echo "$dt $dh" >> "$datafile"
      printf '    n=%-3d dt=%-9s <|dH|>=%-13s acc(meas)=%s\n' "$n" "$dt" "$dh" "$acc"
    done

    npts=$(wc -l < "$datafile")
    if [ "$npts" -lt 3 ]; then
      RESULTS="$RESULTS$cell|FAIL|$expected|n/a|$note
"
      continue
    fi
    pfit=$(awk '{x=log($1); v=$2+0; if (v<=0) next; y=log(v);
                 n++; sx+=x; sy+=y; sxx+=x*x; sxy+=x*y}
                END {if (n>=3) printf "%.3f", (n*sxy-sx*sy)/(n*sxx-sx*sx); else printf "n/a"}' "$datafile")
    verdict="PASS"
    awk -v p="$pfit" -v e="$expected" 'BEGIN{d=p-e; if (d<0) d=-d; exit !(d>0.5)}' && verdict="CHECK"
    [ $fail -eq 1 ] && verdict="PARTIAL"
    printf '    p_fit = %s (expected %s) -> %s\n' "$pfit" "$expected" "$verdict"
    RESULTS="$RESULTS$cell|$verdict|$expected|$pfit|$note
"
  done
done

{
  echo
  echo "============ EQUILIBRATED dH SCALING (8^4 heatbath start, 5+5+5 protocol) ============"
  printf '%-30s %-9s %-10s %-8s %s\n' "cell" "verdict" "expected p" "p_fit" "note"
  printf '%s\n' "--------------------------------------------------------------------------------"
  printf '%s' "$RESULTS" | while IFS='|' read -r c v e p nt; do
    [ -n "$c" ] && printf '%-30s %-9s %-10s %-8s %s\n' "$c" "$v" "$e" "$p" "$nt"
  done
  echo
  echo "Protocol per dt point: load 100-sweep beta=6 heatbath lattice, 5 forced-accept +"
  echo "5 Metropolis (discarded) + 5 measured. <|dH|> of the 5 measured enters the fit."
  echo "CHECK = |p_fit - expected| > 0.5. Data: logs_scaling_therm/<cell>.dat"
} | tee "$LOG_DIR/summary.txt"

#!/usr/bin/env bash
# run_nested_fgi_theory.sh — nested-FGI two-timescale theory verification on 8^4.
#
# Theory (Shcherbakov et al., arXiv:1512.03812 §2.2): with force-gradient
# information at the outer level only and m fixed inner leapfrog substeps,
#
#     dH(h, m) = a*(h/m)^2 + b*h^4 + higher order.
#
# Protocol: single deterministic trajectory per (n_outer, m) grid point —
# identical 100-sweep heatbath start, identical momentum seed, forced accept —
# so every point shares the same start state and hence the same (a, b).
# The signed dH of trajectory 1 is recorded; fit_nested_theory.py then
#   T1: fits (a, b) globally over the (h, m) grid,
#   T2: fits on m <= 3 only and PREDICTS the m in {6, 12} rows,
#   T3: checks h-slice effective exponents (2 -> 4 crossover with m),
#   T4: checks m-slice 1/m^2 decay onto the m-independent b*h^4 plateau.
#
# Run from build/tests:  bash run_nested_fgi_theory.sh
set -u

BIN=./hmc_test
GAUGE="$PWD/therm_8x8_b6.0_hb100.lime"
RES_DIR="$PWD/tunecache_nested_theory"
LOG_DIR="$PWD/logs_nested_theory"
DIM="--dim 8 8 8 8"
TAU=0.5
SEED=12345
TOL=1e-10                    # keep solver noise far below the smallest dH
N_OUTER_GRID="6 8 12 18 27 40"
M_INNER_GRID="1 2 3 6 12"
CELL_TIMEOUT=1800

PREC_MIX="--prec double --prec-sloppy single --prec-precondition single"
ACTION="--dslash-type wilson --mass 0.5"
MG="--mg-levels 2 --mg-block-size 0 4 4 4 4"

[ -x "$BIN" ] || { echo "ERROR: run from build/tests"; exit 1; }
[ -f "$GAUGE" ] || { echo "ERROR: thermalized gauge $GAUGE not found"; exit 1; }
mkdir -p "$RES_DIR" "$LOG_DIR"
if [ -f "$RES_DIR/tunecache.tsv" ] && [ "$BIN" -nt "$RES_DIR/tunecache.tsv" ]; then
  echo "Wiping stale tunecache."
  rm -f "$RES_DIR"/*.tsv
fi

GRID="$LOG_DIR/grid.dat"
{
  echo "# nested-FGI theory grid: single deterministic trajectory per point"
  echo "# wilson 8^4, tau=$TAU, seed=$SEED, tol=$TOL, start=$(basename "$GAUGE")"
  echo "# n_outer  m_inner  h  dH(signed)"
} > "$GRID"

total=0 ok=0
for m in $M_INNER_GRID; do
  for n in $N_OUTER_GRID; do
    total=$((total + 1))
    log="$LOG_DIR/n${n}_m${m}.log"
    # shellcheck disable=SC2086
    QUDA_RESOURCE_PATH="$RES_DIR" timeout "$CELL_TIMEOUT" mpiexec -n 1 "$BIN" $DIM \
      --load-gauge "$GAUGE" \
      $ACTION $PREC_MIX --tol $TOL \
      --hmc-integrator 3 --hmc-n-steps "$n" --hmc-n-inner-steps "$m" --hmc-tau $TAU \
      --hmc-n-trajectories 1 --hmc-thermalization 1 --hmc-momentum-seed $SEED \
      $MG --gtest_filter=HMC.Production > "$log" 2>&1
    rc=$?
    # signed dH of the single (forced-accept) trajectory
    dh=$(grep -E 'hmcRunQuda: (ACCEPTED|REJECTED)' "$log" | awk 'NR==1 {printf "%+.8e", $5+0}')
    if [ $rc -ne 0 ] || [ -z "$dh" ]; then
      echo "  n=$n m=$m : FAILED (rc=$rc), see $log"
      continue
    fi
    ok=$((ok + 1))
    h=$(awk -v t=$TAU -v n=$n 'BEGIN{printf "%.8f", t/n}')
    echo "$n $m $h $dh" >> "$GRID"
    printf '  [%s] n=%-3d m=%-3d h=%-11s dH=%s\n' "$(date '+%H:%M:%S')" "$n" "$m" "$h" "$dh"
  done
done

echo
echo "Grid complete: $ok/$total points -> $GRID"
python3 "$(dirname "$0")/fit_nested_theory.py" "$GRID" | tee "$LOG_DIR/summary.txt"