Merge pull request #50 from E3SM-Project/develop

Develop

methods/islet/Makefile

@@ -28,6 +28,9 @@ pum_sweep: libislet pum_sweep.o
 cslunstab: cslunstab.o
 	$(CXX) cslunstab.o $(LDFLAGS) $(LINK_LAPACK_BLAS) -fopenmp -o cslunstab
 
+pum_perturb_plot: libislet pum_perturb_plot.o
+	$(CXX) pum_perturb_plot.o $(LDFLAGS) $(LINK_LAPACK_BLAS) -L. -lislet -fopenmp -o pum_perturb_plot
+
 clean:
 	rm -f *.o *.so search np4 pum_sweep pum_perturb_plot run_meam1_sweep
 

methods/islet/figures/figs.tex

@@ -129,24 +129,18 @@
 
 % hy figs-methods.hy write-slmmir-script
 % bash run-slmmir-on-basis-lines.sh > slmmir-on-basis-lines-2.txt
-% hy figs-methods.hy plot-slmmir-vs-heuristic # uses slmmir-on-basis-lines-2.txt
-\begin{figure}[tbh]
+% hy figs-methods.hy plot-slmmir-vs-heuristic-ab # uses slmmir-on-basis-lines-2.txt
+\begin{figure}[tb]
   \centering
-  \includegraphics[width=0.48\linewidth]{slmmir-vs-heuristic-gau-nopp-l2}
-  \caption{$l_2$ norm on the nondivergent flow problem
+  \includegraphics[width=1\linewidth]{slmmir-vs-heuristic-ab}
+  \caption{$l_2$ errors for the nondivergent flow problem
     using basis $\basisns_{\np}$ vs.~$a_2(\basisns_{\np})$,
     for a large number of \abtps~bases and $\np=6$ to $10$.
-    The legend lists the marker type for each $\np$.
-    Large red circles outline the bases in Table \ref{tbl:gll}.
-    The configuration uses the Gaussian hills IC and no property preservation.}
-  \label{fig:slmmir-vs-heuristic-a}
-\end{figure}
-\begin{figure}[tbh]
-  \centering
-  \includegraphics[width=0.48\linewidth]{slmmir-vs-heuristic-cos-pp-l2}
-  \caption{Same as Fig.~\ref{fig:slmmir-vs-heuristic-a} except that the configuration
-    uses the cosine bells IC with property preservation.}
-  \label{fig:slmmir-vs-heuristic-b}
+    The legends list the marker type for each $\np$.
+    Large red circles outline the points corresponding to the bases in Table \ref{tbl:gll}.
+    (a) With the Gaussian hills IC and no property preservation.
+    (b) With the cosine bells IC and property preservation.}
+  \label{fig:slmmir-vs-heuristic}
 \end{figure}
 
 % hy figs-adv-diag.hy figs-acc acc-0.txt
@@ -374,3 +368,74 @@
   }
   \label{fig:summit-perf}
 \end{figure}
+
+% hy figs-methods.hy cubed-sphere-subelem-grid-schematic
+\begin{figure}[tb]
+  \centering
+  \includegraphics[width=0.5\linewidth]{cubed-sphere-subelem-grid-schematic}
+  \caption{
+    Cubed-sphere grid (black lines in foreground, gray lines in background)
+    with $\neface\!\times\!\neface$ spectral elements per cube face;
+    $\neface=2$ in this example.
+    The green dashed line outlines the sphere's projection onto the two-dimensional plane of the figure.
+    The upper-right element of the front cubed-sphere face shows the subelement tensor-product Gauss--Lobatto--Legendre (GLL) grid points.
+    In this example, the dynamics solver's subelement grid uses $\npv=4$ (large blue circles) GLL points per dimension,
+    and the transport solver's subelement grid uses $\npt=6$ (small red circles).
+  }
+  \label{fig:cubed-sphere-subelem-grid-schematic}
+\end{figure}
+
+% hy figs-methods.hy isl-1d-schematic
+\begin{figure}[tb]
+  \centering
+  \includegraphics[width=1\linewidth]{isl-1d-schematic}
+  \caption{
+    Illustration of the classical and element-based interpolation semi-Lagrangian methods.
+    See the discussion in Sect.~\ref{sec:setting:sl}.
+  }
+  \label{fig:isl-1d-schematic}
+\end{figure}
+
+% hy figs-methods.hy matrix-schematic
+\begin{figure}[tb]
+  \centering
+  \includegraphics[width=0.5\linewidth]{matrix-schematic}
+  \caption{
+    Correspondence between the ISL space--time operator $\mat{A}$ (top)
+    and the target and source one-dimensional grids (bottom)
+    for one time step of the test problem.
+    In the matrix $\mat{A}$, of which the upper-left corner is pictured,
+    numbered columns correspond to source-grid degrees of freedom (DOF),
+    and numbered rows correspond to target-grid DOF.
+    Nonzeros occur in the red rectangular blocks $\mat{B} \equiv (\mat{\bar{B}} \ \vec{b})$.
+    In this example, the target grid (green dashed line) is advected backward in time
+    so that one subelement grid point moves one element to the left;
+    in the resulting matrix $\mat{A}$, the blocks are shifted one row down.
+  }
+  \label{fig:matrix-schematic}
+\end{figure}
+
+% bash run-instab-imgs.sh
+% hy figs-adv-diag.hy fig-instab
+\begin{figure}[tb]
+  \centering
+  \includegraphics[width=0.75\linewidth]{img-instab}
+  \caption{
+    Images for unstable (top) and stable (bottom) ISL transport.
+    The problem is nondivergent flow with the slotted cylinders IC,
+    with $\neface = 20$, $\npv = \npt = 6$, and the large time step.
+    The snapshot is at the end of day 11 of the first cycle;
+    the images are zoomed to just the region containing the slotted cylinders.
+    The color range is [-0.05, 1.15],
+    which clips the top-left image's range of [-30.4, 32.7].
+    The bases are GLL (top) and Islet (bottom),
+    without (left) and with (right) property preservation.
+    The GLL basis yields an unstable ISL method.
+    Although a nonlinear property preservation step makes the method stable,
+    the linear advection operator's instability still manifests as spurious oscillations.
+    The Islet basis yields a stabilized ISL method;
+    the nonlinear property preservation step now just controls mass conservation and extrema,
+    as intended.
+  }
+  \label{fig:islet-vs-gll-img}
+\end{figure}

methods/islet/figures/run-instab-imgs.sh

@@ -0,0 +1,22 @@
+exe=../../slmm/slmmir
+
+ne=20
+np=6
+nstep=$(($ne * 6))
+we=110
+
+glbcdrs=(none caas-node none caas-node)
+bases=(Gll Gll GllNodal GllNodal)
+ncycles=(1 1 1 1)
+
+for i in $(seq 0 3); do
+    ncycle=${ncycles[$i]}
+    glbcdr=${glbcdrs[$i]}
+    basis=${bases[$i]}
+    echo $i $ncycle $glbcdr $basis
+    name=instab-nondiv-ne${ne}np${np}-${glbcdr}-cyc${ncycle}-$basis
+    time=$(($ncycle * 12))
+    cmd="OMP_PROC_BIND=false OMP_NUM_THREADS=8 $exe -method pcsl -ode nondivergent -ic gaussianhills -ic cosinebells -ic slottedcylinders -T $time -nsteps $nstep -timeint exact -ne $ne -np $np -dmc eh -d2c -mono $glbcdr -lim caas -we $we -io-type internal -io-start-cycle $ncycle -res 256 -o $name -rit -prefine 0 -basis $basis -rhot0"
+    echo "cmd> $cmd"
+    eval "$cmd"
+done

methods/islet/pum_perturb_plot.cpp

@@ -0,0 +1,51 @@
+#include "islet_isl.hpp"
+#include "islet_pum.hpp"
+
+#include <cstdio>
+
+static void write_array (const std::string& name, const std::vector<Real>& a) {
+  printf("%s [", name.c_str());
+  for (const auto& e : a) printf(" %1.5e", e);
+  printf("]\n");
+}
+
+static void gen_plot_data (const islet::Operator::ConstPtr& op, const Int np,
+                           const bool print_perturb = true) {
+  static const Int nperturb = 48;
+  std::vector<Real> perturb(nperturb), meam1(nperturb);
+  for (Int i = 0; i < nperturb; ++i) perturb[i] = 0.1/std::pow(1.18, i);
+  if (print_perturb) write_array(">> perturb", perturb);
+  const auto oim = std::make_shared<islet::OperatorInterpMethod>(np, op);
+  const auto im = std::make_shared<InterpMethod>(oim);
+# pragma omp parallel for
+  for (Int i = 0; i < nperturb; ++i) {
+    meam1[i] = 0;
+    pum::Options o;
+    o.threaded = false;
+    o.ntrial = 33;    
+    o.perturb = perturb[i];
+    for (const Int ne : {4, 7, 15})
+      for (const Int mec_ne : {33}) {
+        o.ne = ne;
+        o.mec_ne = mec_ne;
+        pum::PerturbedUniformMeshMetric pum(im, o);
+        meam1[i] = std::max(meam1[i], pum.run());
+      }
+  }
+  write_array(">> meam1", meam1);
+}
+
+int main (int argc, char** argv) {
+  const auto gll_best = islet::Operator::create(islet::Operator::gll_best);
+  const auto uofs = islet::Operator::create(islet::Operator::uniform_offset_nodal_subset);
+  bool first = true;
+  for (int np = 2; np <= islet::np_max; ++np) {
+    if (np >= 4) {
+      printf(">>> gll_best %d\n", np);
+      gen_plot_data(gll_best, np, false);
+    }
+    printf(">>> uniform_offset_nodal_subset %d\n", np);
+    gen_plot_data(uofs, np, first);
+    first = false;
+  }
+}

methods/islet/readme.txt

@@ -53,15 +53,15 @@ installation, then
     ln -s make.inc.gnu make.inc
     make -j16
 Optionally run regression tests:
-    OMP_NUM_THREADS=16 KMP_AFFINITY=balanced python2 slmm_runtests.py
+    OMP_NUM_THREADS=16 KMP_AFFINITY=balanced python slmm_runtests.py
 Bash scripts in the methods/islet/figures directory call the slmmir program.
 
 We use the language hy to create the figures. hy is a Lisp that compiles to
 Python AST. We used hy 0.18.0 ('pip install hy' for the latest version) with
 CPython 3.7.6 provided by Anaconda 3.
 
-The code used to obtain performance data on Summit will be part of main E3SM
-soon. The exact version used to generate the data is archived here:
+The code used to obtain performance data on Summit is part of the main E3SM
+repo. The exact version used to generate the data is archived here:
     https://github.com/ambrad/E3SM/releases/tag/islet-2d-paper-summit-sl-gpu-timings
 The data are here:
     https://github.com/E3SM-Project/perf-data/tree/main/nhxx-sl-summit-mar2021