Facilitate the making of the incompressible convergence plot (#147)

The main addition is pyro/incompressible/tests/convergence_error.py which automatically generates the file which contains the errors at each resolution. Perhaps some cleanup/syntax check would be good.

Author: simonguichandut

docs/source/incompressible_basics.rst

@@ -62,7 +62,18 @@ shown below are run as:
    pyro_sim.py incompressible converge inputs.converge.128 vis.dovis=0
 
 The error is measured by comparing with the analytic solution using
-the routine ``incomp_converge_error.py`` in ``analysis/``.
+the routine ``incomp_converge_error.py`` in ``analysis/``. To generate
+the plot below, run
+
+.. prompt:: bash
+
+   python incompressible/tests/convergence_errors.py convergence_errors.txt
+
+or ``convergence_errors_no_limiter.txt`` after running with that option. Then:
+
+.. prompt:: bash
+
+   python incompressible/tests/convergence_plot.py
 
 .. image:: incomp_converge.png
    :align: center

pyro/analysis/incomp_converge_error.py

@@ -16,36 +16,42 @@
 """
 
 
-if len(sys.argv) != 2:
-    print(usage)
-    sys.exit(2)
+def get_errors(file):
 
+    sim = io.read(file)
+    myd = sim.cc_data
+    myg = myd.grid
 
-try:
-    file1 = sys.argv[1]
-except IndexError:
-    print(usage)
-    sys.exit(2)
+    # numerical solution
+    u = myd.get_var("x-velocity")
+    v = myd.get_var("y-velocity")
 
-sim = io.read(file1)
-myd = sim.cc_data
-myg = myd.grid
+    t = myd.t
 
-# numerical solution
-u = myd.get_var("x-velocity")
-v = myd.get_var("y-velocity")
+    # analytic solution
+    u_exact = myg.scratch_array()
+    u_exact[:, :] = 1.0 - 2.0*np.cos(2.0*math.pi*(myg.x2d-t))*np.sin(2.0*math.pi*(myg.y2d-t))
 
-t = myd.t
+    v_exact = myg.scratch_array()
+    v_exact[:, :] = 1.0 + 2.0*np.sin(2.0*math.pi*(myg.x2d-t))*np.cos(2.0*math.pi*(myg.y2d-t))
 
-# analytic solution
-u_exact = myg.scratch_array()
-u_exact[:, :] = 1.0 - 2.0*np.cos(2.0*math.pi*(myg.x2d-t))*np.sin(2.0*math.pi*(myg.y2d-t))
+    # error
+    udiff = u_exact - u
+    vdiff = v_exact - v
 
-v_exact = myg.scratch_array()
-v_exact[:, :] = 1.0 + 2.0*np.sin(2.0*math.pi*(myg.x2d-t))*np.cos(2.0*math.pi*(myg.y2d-t))
+    return udiff.norm(), vdiff.norm()
 
-# error
-udiff = u_exact - u
-vdiff = v_exact - v
 
-print("errors: ", udiff.norm(), vdiff.norm())
+if __name__ == "__main__":
+    if len(sys.argv) != 2:
+        print(usage)
+        sys.exit(2)
+
+    try:
+        file = sys.argv[1]
+    except IndexError:
+        print(usage)
+        sys.exit(2)
+
+    errors = get_errors(file)
+    print("errors: ", errors[0], errors[1])

pyro/incompressible/tests/convergence_errors.py

@@ -0,0 +1,37 @@
+# This should be run from /pyro (where the output files are located)
+import glob
+import sys
+
+from pyro.analysis.incomp_converge_error import get_errors
+
+
+def create_file(filename="convergence_errors.txt"):
+
+    # Parse all converge*.h5 outputs
+    # store the last file at each resolution
+    fnames = glob.glob("converge*.h5")
+    res = set([f.split("_")[1] for f in fnames])
+    res = list(res)
+    res.sort(key=int)
+
+    simfiles = []
+    for r in res:
+        fnames = glob.glob(f"converge_{r}_*.h5")
+        last = max([int(f.split("_")[-1].split(".")[0]) for f in fnames])
+        simfiles.append(f"converge_{r}_{last:04d}.h5")
+
+    # Create the file
+    with open("incompressible/tests/" + filename, "w") as f:
+        f.write("# convergence of incompressible converge problem\n")
+        f.write("# (convergence measured with analysis/incomp_converge_error.py)\n")
+        f.write("#\n")
+        for r, file in zip(res, simfiles):
+            errors = get_errors(file)
+            f.write(f"{r.rjust(3)} {errors[0]:.14f} {errors[1]:.14f}\n")
+
+
+if __name__ == "__main__":
+    if len(sys.argv) == 2:
+        create_file(sys.argv[1])
+    else:
+        create_file()

…incompressible/tests/incomp_converge.txt …ompressible/tests/convergence_errors.txtpyro/incompressible/tests/incomp_converge.txt renamed to pyro/incompressible/tests/convergence_errors.txt pyro/incompressible/tests/incomp_converge.txt renamed to pyro/incompressible/tests/convergence_errors.txt

@@ -1,9 +1,6 @@
 # convergence of incompressible converge problem
 # (convergence measured with analysis/incomp_converge_error.py)
 #
- 32 0.0220739241002 0.0220739241002
+ 32 0.02207392410018 0.02207392410018
  64 0.00682161894998 0.00682161894998
-128 0.00214105266633 0.00214105266633
-
-
-
+128 0.00214081341364 0.00214081341364

pyro/incompressible/tests/convergence_errors_no_limiter.txt

@@ -0,0 +1,6 @@
+# convergence of incompressible converge problem
+# (convergence measured with analysis/incomp_converge_error.py)
+#
+ 32 0.01491254251261 0.01491254251261
+ 64 0.00216262704768 0.00216262704768
+128 0.00036281268740 0.00036281268740

pyro/incompressible/tests/convergence_plot.py

@@ -1,20 +1,37 @@
+import os
+
 import matplotlib.pyplot as plt
 import numpy as np
 
 
 def plot_convergence():
 
-    data = np.loadtxt("incomp_converge.txt")
-
-    nx = data[:, 0]
-    errlim = data[:, 1]
-
     ax = plt.subplot(111)
     ax.set_xscale('log')
     ax.set_yscale('log')
 
-    plt.scatter(nx, errlim, marker="x", color="r", label="limiting")
-    plt.plot(nx, errlim[0]*(nx[0]/nx)**2, "--", color="k")
+    files = ["incompressible/tests/convergence_errors.txt",
+             "incompressible/tests/convergence_errors_no_limiter.txt"]
+
+    if not os.path.exists(files[0]):
+        print("Could not find convergence_errors.txt")
+    else:
+        data = np.loadtxt(files[0])
+        nx = data[:, 0]
+        errlim = data[:, 1]
+
+        plt.scatter(nx, errlim, marker="x", color="r", label="limiting")
+        plt.plot(nx, errlim[0]*(nx[0]/nx)**2, "--", color="k")
+
+    if not os.path.exists(files[1]):
+        print("Could not find convergence_errors_no_limiter.txt")
+    else:
+        data = np.loadtxt(files[1])
+        nx = data[:, 0]
+        errlim = data[:, 1]
+
+        plt.scatter(nx, errlim, marker="x", color="b", label="limiting")
+        plt.plot(nx, errlim[0]*(nx[0]/nx)**2, "--", color="k")
 
     plt.xlabel("number of zones")
     plt.ylabel("L2 norm of abs error")