Merge pull request #51 from ambrad/ambrad/islet2-methods-update

Islet2 methods update

Author: ambrad

methods/islet/figures/figs-adv-diag.hy

@@ -440,7 +440,8 @@
 (defn norm->str [n]
   (get {:l1 "$l_1$" :l2 "$l_2$" :li "$l_{\infty}$"} n))
 
-(defn plot1-slmmir-vs-heuristic [c d nps prop-preserve ic norm pum-thr lebesgue]
+(defn plot1-slmmir-vs-heuristic [c d nps prop-preserve ic norm pum-thr lebesgue
+                                 &optional [jcp False]]
   (sv npa npy.array fs 11
       plot (if lebesgue pl.semilogy pl.loglog))
   (for [np nps]
@@ -466,7 +467,10 @@
                    :fontsize fs)))
   (my-grid)
   (pl.title (.format (+ "Test problem vs.~heuristic:\n"
-                        "nondivergent flow, 1.5$^\circ$, {}, long steps\n"
+                        "nondivergent flow, "
+                        (if jcp "$n_e = 20$, " "1.5$^\circ$, ")
+                        "{}, "
+                        (if jcp "120 steps\n" "long steps\n")
                         "$p$-refinement, {}")
                      (futils.ic-short2long ic)
                      (+ (if prop-preserve "" "no ") "property preservation"))
@@ -498,7 +502,7 @@
       (sv t (nth norm-pp-ic-tuples i))
       (pl.subplot 1 2 (inc i))
       (plot1-slmmir-vs-heuristic c d nps (nth t 1) (nth t 2) (nth t 0) pum-thr
-                                 lebesgue)
+                                 lebesgue :jcp True)
       (sv f (pl.gcf))
       (f.text (nth (, 0.02 0.52) i) 0.05 (nth (, "(a)" "(b)") i)
               :fontsize fs))))
@@ -644,7 +648,7 @@
 ;;; drivers
 
 (when-inp ["dev-parse"]
-  (sv fname "data/search-0.txt"
+  (sv fname "data/mar21/search-0.txt"
       blns (parse-search-list fname))
   (print (len blns))
   (for [b blns]
@@ -664,7 +668,7 @@
                      "search-findnodal_given_bestosn-7.txt")
       blns [])
   (for [fname data-fnames]
-    (.extend blns (parse-search-list (+ "data/" fname))))
+    (.extend blns (parse-search-list (+ "data/mar21/" fname))))
   (sv blns (uniquify-search-list blns))
   (write-slmmir-script blns script-fname))
 
@@ -674,7 +678,7 @@
       lebesgue False
       d {})
   (for [fname fnames]
-    (sv d (parse-slmmir-output (+ "data/" fname) :d d :lebesgue lebesgue)))
+    (sv d (parse-slmmir-output (+ "data/mar21/" fname) :d d :lebesgue lebesgue)))
   (for [(, norm pp ic) (, (, :l2 False "gau") (, :l2 True "cos"))]
     (plot-slmmir-vs-heuristic
       c d (.format (+ "{}slmmir-vs-heuristic-{}-{}-{}" (if lebesgue "-leb" ""))
@@ -688,7 +692,7 @@
       lebesgue False
       d {})
   (for [fname fnames]
-    (sv d (parse-slmmir-output (+ "data/" fname) :d d :lebesgue lebesgue)))
+    (sv d (parse-slmmir-output (+ "data/mar21/" fname) :d d :lebesgue lebesgue)))
   (plot-slmmir-vs-heuristic-ab
    c d (+ c.fig-dir "slmmir-vs-heuristic-ab")
    (, (, :l2 False "gau") (, :l2 True "cos"))
@@ -702,13 +706,13 @@
 
 (when-inp ["pum-vs-perturb"]
   (sv fname "pum_perturb_plot-041021.txt"
-      d (parse-pum-vs-perturb (+ "data/" fname))
+      d (parse-pum-vs-perturb (+ "data/mar21/" fname))
       c (futils.get-context))
   (plot-pum-vs-perturb c d (+ c.fig-dir "pum-vs-perturb")))
 
 (when-inp ["meam1-and-pum-vs-dx"]
   (sv c (futils.get-context)
-      data-dir "data/"
+      data-dir "data/mar21/"
       meam1-fname (+ data-dir "run_meam1_sweep-np8.txt")
       method-fnames (zip (, "gll_best" "gll_natural" "uniform_offset_nodal_subset")
                          (lfor fname (, "pum_sweep-np8-gll_best.txt"

methods/islet/figures/figs-methods.hy

@@ -0,0 +1,120 @@
+(require [amb3 [*]])
+(import [amb3 [*]]
+        [figsutils [*]]
+        math colorsys)
+
+(assoc matplotlib.rcParams "savefig.dpi" 300)
+(do (pl-require-type1-fonts))
+
+(defn get-errs [d ne np cycle]
+  (get d "exact" "movingvortices" 1 "pcsl" "none" "none" 0 np "vor" cycle ne :C))
+
+(defn make-vortex-cmap []
+  (sv lsc (fn [colors lo hi n]
+            ((matplotlib.colors.LinearSegmentedColormap.from-list
+              "filament" colors 11)
+             (npy.linspace lo hi n)))
+      h2r colorsys.hsv-to-rgb
+      s 0.2
+      red (lsc [(h2r 0 1 1) (h2r 0 s 1)] 0 1 11)
+      b (/ 2 3)
+      blue (lsc [(h2r b 1 1) (h2r b s 1)] 1 0 11)
+      cmap (matplotlib.colors.ListedColormap
+            (npy.concatenate (, blue red) :axis 0)))
+  cmap)
+
+(defn fig-vortex [c d configs]
+  (print "TODO config labels")
+  (defn draw-circle []
+    (sv th (npy.linspace 0 (* 2 math.pi) 512)
+        r 0.995
+        x (* r (npy.cos th)) y (* r (npy.sin th))
+        g 0.8)
+    (pl.plot x y "-" :lw 0.75 :color [g g g]))
+  (sv ncol (len configs)
+      fs 11
+      cmap (make-vortex-cmap)
+      ecmap (matplotlib.colors.LinearSegmentedColormap.from-list
+             "filament" [[1 1 1] [0 0 0]]))
+  (with [(pl-plot (, (* 2 ncol) 4) (+ c.fig-dir "vortex")
+                  :format "png"
+                  :tight False)]
+    (for [(, icon config) (enumerate configs)]
+      (sv (, ne np nstep cycle recon) config
+          e (get-errs d ne np cycle)
+          arrs (read-slmmir-io-arrays
+                (+ c.data-dir f"vortex-imgs/ne{ne}np{np}nstep{nstep}{recon}.bin")))
+      (sv img (get arrs (+ 1 (* 3 cycle)))
+          diff (/ (npy.abs (- img (get arrs (+ 2 (* 3 cycle)))))
+                  (npy.max (npy.abs (get arrs (+ 2 (* 3 cycle))))))
+          imgs [img diff])
+      (for [i (range 2)]
+        (sv img (get imgs i)
+            mask (!= img 0)
+            ext (, (npy.min (get img mask)) (npy.max img)))
+        (when (zero? i)
+          (sv nmask (npy.logical-not mask)))
+        (sv (get img nmask) npy.nan)
+        (print "corner" (get img (, 0 0)))
+        (print "ext" ext)
+        (if (zero? i)
+          (sv vmin 0.45 vmax 1.55)
+          (sv vmax (get e :li)
+              vmin 0))
+        (sv rect [(/ icon ncol) (- 1 (* i 0.5))
+                  (/ 0.97 ncol) (* 0.97 0.5)]
+            ax (pl.axes rect))
+        (pl.imshow img :cmap (if (zero? i) cmap ecmap)
+                   :vmin vmin :vmax vmax
+                   :origin "upper"
+                   :extent [-1 1 -1 1])
+        (pl.axis "off")
+        (pl.text -1 -1
+                 (.format "({})" (if (zero? i)
+                                   (get "abcd" icon)
+                                   (get "efgh" icon)))
+                 :ha "left" :va "bottom" :fontsize fs)
+        (draw-circle)
+        (unless (zero? i)
+          (pl.text 0 0.5
+                   (.format "$n_e \, {}$, $n_p \, {}$, cycle {}, $n_{{step}} \, {}$\n"
+                            ne np cycle (* cycle nstep))
+                   :ha "center" :fontsize fs)
+          (unless (= recon "constant")
+            (for [k (range 3)]
+              (pl.text -0.4 (+ 0.0 (* 0.2 (- 2 k)))
+                       (.format (+ (get ["$l_1$" "$l_2$" "$l_{{\infty}}$"] k)
+                                   " {:1.2e}")
+                                (get e (get (, :l1 :l2 :li) k)))
+                       :ha "right" :fontsize fs)))
+          (when (= recon "constant")
+            (pl.text 0 0.9 "constant in subcell" :ha "center" :fontsize fs)))
+        (cond [(and (zero? i) (= icon 0))
+               (sv f (pl.gcf)
+                   cax (f.add-axes [(- (get rect 0) 0.04)
+                                    (+ (get rect 1) (* 0.1 (get rect 3)))
+                                    (* 0.1 (get rect 2))
+                                    (* 0.8 (get rect 3))]))
+               (pl.colorbar :cax cax :ticks [0.5 0.75 1 1.25 1.5]
+                            :aspect 10)
+               (cax.yaxis.set-ticks-position "left")]
+              [(= i 1)
+               (sv f (pl.gcf)
+                   cax (f.add-axes [(+ (get rect 0) (* 0.1 (get rect 2)))
+                                    (+ (get rect 1) 0.1)
+                                    (* 0.8 (get rect 2))
+                                    (* 0.05 (get rect 3))])
+                   cb (pl.colorbar :cax cax :orientation "horizontal"
+                                   :ticks [0 vmax]
+                                   :shrink 0.6 :aspect 10))
+               (cb.ax.set-xticklabels ["0" (.format "{:1.2e}" vmax)])])))))
+
+(when-inp ["fig-vortex"]
+  (sv c (jcp-context (get-context))
+      d (acc-parse (+ c.data-dir "vortex-imgs.txt")
+                   :map-nstepfac jcp-nenp2nstepfac)
+      configs [(, 3 13 72 1 "constant")
+               (, 3 13 72 1 "bilin")
+               (, 6 13 144 1 "bilin")
+               (, 12 13 288 2 "bilin")])
+  (fig-vortex c d configs))

methods/islet/figures/figs-vortex.hy

@@ -3,8 +3,8 @@
 
 (defn get-context []
   (sv c (Box)
-      c.data-dir "data/"
-      c.fig-dir "figs/"
+      c.data-dir "data/mar21/"
+      c.fig-dir "figs/drafts/"
       c.odes (, "rotate" "divergent" "nondivergent")
       c.cdrs (, (, "none" "none") (, "caas-node" "caas") (, "caas" "caas"))
       c.nstepfacs (, 1 5)
@@ -14,32 +14,88 @@
       c.nes (, 5 10 20 40 80)
       c.nps (, 4 6 8 9 12) ;(, 4 5 6 7 8 9 10 11 12 13 16)
       c.ics (, "gau" "cos" "slo")
-      c.npclrs  {4 "g" 5 "m" 6 "r" 7 "c" 8 "k" 9 "b" 10 "g" 11 "c" 12 "m" 13 "m" 16 "g"}
+      c.npclrs  {4 "g" 5 "m" 6 "r" 7 "c" 8 "k" 9 "b" 10 "g" 11 "c" 12 "r" 13 "m" 16 "g"}
       c.npmarks {4 "o" 5 "x" 6 "s" 7 "x" 8 "p" 9 "+" 10 "." 11 "^" 12 "." 13 "*" 16 "."})
   c)
 
 (defn flow-short2long [flow]
   (get {"divergent" "divergent flow"
         "nondivergent" "nondivergent flow"
-        "rotate" "solid-body rotation"} flow))
+        "rotate" "solid-body rotation"
+        "movingvortices" "moving vortices"} flow))
 
 (defn ic-short2long [ic]
   (get {"gau" "Gaussian hills"
         "cos" "cosine bells"
         "slo" "slotted cylinders"} ic))
 
-(defn nes->degstrs [nes]
+(defn nes->degstrs [nes &optional [convert-all False]]
   (sv x [] xstr [])
   (for [ne nes]
-    (sv deg (geton {5 6 10 3 20 "1.5" 40 "0.75" 80 "0.375" 160 "0.1875"} ne))
-    (when (none? deg) (continue))
+    (sv deg (geton {2 15 5 6 10 3 20 "1.5" 40 "0.75" 80 "0.375" 160 "0.1875"} ne))
     (.append x ne)
-    (.append xstr (.format "${}^{{\circ}}$" deg)))
+    (if (none? deg)
+      (do (sv deg (/ 30 ne))
+          (.append xstr (.format "${:1.4f}^{{\circ}}$" deg)))
+      (.append xstr (.format "${}^{{\circ}}$" deg))))
   {:ne x :deg xstr})
 
 (defn cdr-name [short]
   (get {"caas" "CAAS-CAAS" "caas-node" "CAAS-point"} short))
 
+;;; parse cmd, L, M, C slmmir output
+
+(defn acc-parse [fname &optional map-nstepfac]
+  (sv txt (.split (readall fname) "\n")
+      bo None d {})
+  (for [ln txt]
+    (sv ln2 (cut ln 0 2))
+    (cond [(and (= "cm" ln2) (in "cmd>" ln))
+           (sv cmd ln c (parse-cmd ln :map-nstepfac map-nstepfac))]
+          [(= ln2 "M ")
+           (sv mp (parse-midpoint-check ln))
+           (unless (none? (geton d #*(+ (cmd->key-base c)
+                                        (, (:ic mp) cyc (:ne c) :M))))
+             ;; handle repeated ic used for src-term OOA testing.
+             (assoc mp :ic (+ (:ic mp) "-src")))
+           (assoc-nested d (+ (cmd->key-base c) (, (:ic mp) cyc (:ne c) :M))
+                         {:l1 (:l1 mp) :l2 (:l2 mp)})]
+          [(= ln2 "C ")
+           (cond [(in "cycle" ln) (sv (, - - cyc) (sscanf ln "s,s,i"))]
+                 [(or (in "PASS" ln) (in "FAIL" ln))
+                  (sv (, - pf) (sscanf ln "s,s"))
+                  (when (and (= pf "FAIL") (!= (:mono c) "none") (<= cyc 10))
+                    (prf "FAIL {}" cmd))]
+                 [:else
+                  (sv cl (parse-C ln))
+                  (unless (none? (geton d #*(+ (cmd->key-base c)
+                                               (, (:ic cl) cyc (:ne c) :C))))
+                    ;; handle repeated ic used for src-term OOA testing.
+                    (assoc cl :ic (+ (:ic cl) "-src")))
+                  #_(print "found:" (+ (cmd->key-base c) (, (:ic cl) cyc (:ne c) :C)))
+                  (assoc-nested d (+ (cmd->key-base c) (, (:ic cl) cyc (:ne c) :C))
+                                cl)])]
+          [(= ln2 "L ")
+           (cond [(in "L file" ln)
+                  (assoc-nested d (+ (cmd->key-base c)
+                                     (, "cos" cyc (:ne c) :L :mixing-file))
+                                (last (.split ln)))
+                  ;; If nothing else has marked bo as done, do it now.
+                  (unless (none? bo)
+                    (assoc bo :done True))]
+                 [(in "phimin" ln)
+                  (sv (, - ic - l1 - l2 - li - phimin - phimax)
+                      (sscanf ln "s,s,s,f,s,f,s,f,s,f,s,f"))
+                  (assoc-nested d (+ (cmd->key-base c) (, ic cyc (:ne c) :Lerr))
+                                {:l1 l1 :l2 l2 :li li :phimin phimin :phimax phimax})]
+                 [:else
+                  (sv bo (parse-bakeoff-diag bo ln (:timeint c)))])
+           (when (and (not (none? bo)) (:done bo))
+             (for [(, k v) (.items bo)]
+               (assoc-nested d (+ (cmd->key-base c) (, "cos" cyc (:ne c) :L k)) v))
+             (sv bo None))]))
+  d)
+
 ;;; slmmir I/O
 
 (defn read-slmmir-io-arrays [fname &optional beg end stride]
@@ -91,6 +147,7 @@
 
 ;;; parse slmmir text output
 
+;; (map-nstepfac nstepfac-initial ne np)
 (defn parse-cmd [cmd &optional map-nstepfac]
   (sv toks (.split cmd))
   (defn int-or-none [x]
@@ -105,7 +162,8 @@
   (for [e (.items keys)]
     (assoc d (keyword (first e)) ((second e) (get-key-val (first e)))))
   (sv nstepfac (/ (:nsteps d) (:ne d) 6))
-  (unless (none? map-nstepfac) (sv nstepfac (map-nstepfac nstepfac)))
+  (unless (none? map-nstepfac)
+    (sv nstepfac (map-nstepfac nstepfac (:ne d) (:np d))))
   (assoc d :nstepfac (int nstepfac))
   (when (= (:timeint d) "exact") (assoc d :prefine 0))
   d)
@@ -224,3 +282,32 @@
   (if (in " offst " ln)
       (parse-search-offset-nodal-subset-line ln)
       (parse-search-nodal-subset-line ln)))
+
+;;; JCP manuscript revision utilities
+
+(defn jcp-context [c]
+  (sv c.data-dir "data/feb24/"
+      c.nps (, 4 6 8 11 13))
+  c)
+
+;; Compute ne for np given tne for np=4 such that (ne,np) has the maximal #DOF
+;; <= $DOF for (tne,4).
+(defn jcp-tne2ne [tne np]
+  (// (* tne 3) (dec np)))
+
+;; The floating point actual value of resolution converted to np=4.
+(defn jcp-effective-tne [tne np]
+  (/ (* (jcp-tne2ne tne np) (dec np)) 3))
+
+(defn jcp-np4-ne [ne np]
+  (/ (* ne (dec np)) 3))
+
+;; Deduce nstepfac 1 or 5 from parse-cmd's internal values.
+(defn jcp-nenp2nstepfac [nstepfac ne np]
+  (sv nsteps (* nstepfac ne 6)
+      tne (/ (* ne (dec np)) 3)
+      tnstepfac (/ nsteps tne 6))
+  (cond [(< tnstepfac 3) 1]
+        [(< tnstepfac 8) 5]
+        [:else (raisefmt "jcp-nenp2nstepfac failed: {} {} {} => {}"
+                         nstepfac ne np tnstepfac)]))

methods/islet/figures/figsutils.hy

@@ -0,0 +1,28 @@
+cat $0
+
+exe=slmmir
+
+ctr=0
+function run {
+    ctr=$(expr $ctr + 1)
+    cmd="OMP_NUM_THREADS=64 KMP_AFFINITY=balanced ../$exe -method pcsl -ic gaussianhills -ic cosinebells -ic correlatedcosinebells -ic slottedcylinders -we 0 -rit -T 12 -d2c -ode $ode -ne $ne -np $np -nsteps $nstep -prefine 0 -mono none -lim none -timeint exact -rotate-grid"
+    grepcmd='grep "^C \|^L \|^M "'
+    echo "cmd> $ctr $cmd"
+    eval "$cmd | $grepcmd"
+}
+
+for tne in 5 10 20 40 80 160; do
+    for nstepfac in 1 5; do
+        for ode in nondivergent; do
+            tgtres=$(($tne * 3))
+            for np in $(seq 4 13) 16; do
+                ne=$(($tgtres / ($np - 1)))
+                nstep=$(($nstepfac * 2 * ($np - 1) * $ne))
+                if [[ $ne == 1 ]]; then
+                    continue
+                fi
+                run
+            done
+        done
+    done
+done