convergence utility tools (#165)

This pr makes convergence.py more general so that it takes any variable we want with the arbitrary multiplicative factor of resolutions apart between plot files. The default is still density and twice the resolution between the fine plot files and the coarse plot files.

convergence_plot.py allows users to give 3 or more output files which the corresponding resolution differs by a constant multiplicative constant.

The goal here is to duplicate the Richardson utility tool that we had in castro/amrex. It outputs the table that tells us the convergence rate as well as a plot.

Author: zhichen3

docs/source/analysis.rst

@@ -23,10 +23,31 @@ with their data.
 * ``analysis/``
 
   * ``convergence.py``: this compares two files with different
-    resolutions (one a factor of 2 finer than the other).  It coarsens
-    the finer data and then computes the norm of the difference.  This
-    is used to test the convergence of solvers.
-
+    resolutions (one a factor of N-times, (default is 2 times or ``N=2``) finer
+    than the other). It coarsens the finer data and then computes the
+    norm of the difference.  The default variable being compared is
+    density, ``var_name=density`` by default. but it takes arbitary variable name as input.
+    This is used to test the convergence of solvers.
+
+    usage: ``./convergence.py fine_file coarse_file var_name N``
+
+  * ``convergence_plot.py``: this compares 3 or more files with different
+    resolutions (factor of N-times (default is 2 times) finer than the
+    others) by using ``convergence.py``. It prints out a table that summarizes
+    the L-2 Norm of all the variables by coarsening the finer data and comparing
+    it to the coarse data, as well as the order of convergence between the
+    L-2 Norms. It also produces a plot in the end to compare the theoretical
+    order of convergence to the actual order of convergence from the datasets.
+    The data files should enter from the finest file to the coarsest file. There
+    should be at least 3 files, ``fine_file med_file coarse_file``.
+    The output file name for the convergence plot, ``out = convergence_plot.pdf``,
+    by default. The theoretical order of convergence for the solver, ``order=2``,
+    by default. The multiplicative factor between each data files are, ``resolution=2``,
+    by default. This is used to test the order of the convergence of the solver
+    when there is no analytical solutions.
+
+    usage: ``./convergence_plot.py files* -o out -n order -r resolution``
+  
   * ``dam_compare.py``: this takes an output file from the
     shallow water dam break problem and plots a slice through the domain
     together with the analytic solution (calculated in the script).

pyro/analysis/convergence.py

@@ -6,39 +6,47 @@
 
 import pyro.util.io_pyro as io
 
-# read in two files -- one twice the resolution of the other, and
+# read in two files -- one N-times the resolution of the other, and
 # compute the error by averaging down
 
 usage = """
-      usage: ./convergence.py fine coarse
+      usage: ./convergence.py fine coarse variable_name[optional, default=density] N[default=2]
 """
 
 
-def compare(fine, coarse):
+def compare(fine, coarse, var_name, N):
 
-    dens = coarse.get_var("density")
-    dens_avg = fine.restrict("density", N=2)
+    var = coarse.get_var(var_name)
+    var_avg = fine.restrict(var_name, N=N)
 
     e = coarse.grid.scratch_array()
-    e.v()[:, :] = dens.v() - dens_avg.v()
+    e.v()[:, :] = var.v() - var_avg.v()
 
     return float(np.abs(e).max()), e.norm()
 
 
 def main():
-    if len(sys.argv) != 3:
+    if (len(sys.argv) > 5 or len(sys.argv) < 3):
         print(usage)
         sys.exit(2)
 
     fine = sys.argv[1]
     coarse = sys.argv[2]
 
+    var_name = "density"
+    N = 2
+    if len(sys.argv) > 3:
+        var_name = sys.argv[3]
+
+        if len(sys.argv) > 4:
+            N = int(sys.argv[4])
+
     ff = io.read(fine)
     cc = io.read(coarse)
 
-    result = compare(ff.cc_data, cc.cc_data)
+    result = compare(ff.cc_data, cc.cc_data, var_name, N)
 
-    print("inf norm of density: ", result)
+    print(f"inf norm and L2 norm of {var_name}: ", result)
 
 
 if __name__ == "__main__":

pyro/analysis/convergence_plot.py

@@ -0,0 +1,129 @@
+#!/usr/bin/env python3
+
+'''
+This file prints summary of the convergence given at least 3 plot files
+that are differ by a constant multiplicative resolution factor.
+It prints out a table as well as a plot.
+'''
+
+import argparse
+import sys
+
+import convergence
+import matplotlib.pyplot as plt
+import numpy as np
+from prettytable import PrettyTable
+
+import pyro.util.io_pyro as io
+
+DESCRIPTION = "This file does convergence plotting for the solvers"
+INPUT_HELP = '''The location of the input files: enter from finest to coarset files,
+at least 3 files are required'''
+OUTPUT_HELP = "File name for the convergence plot"
+ORDER_HELP = "The theoretical order of convergence for the solver"
+RESOLUTION_HELP = "Multiplicative Resolution between input files, default is 2"
+
+parser = argparse.ArgumentParser(description=DESCRIPTION)
+parser.add_argument("input_file", nargs='+', help=INPUT_HELP)
+parser.add_argument("-o", "--out", default="convergence_plot.pdf", help=OUTPUT_HELP)
+parser.add_argument("-n", "--order", default=2, type=int, help=ORDER_HELP)
+parser.add_argument("-r", "--resolution", default=2, type=int, help=RESOLUTION_HELP)
+
+args = parser.parse_args(sys.argv[1:])
+
+
+def convergence_plot(dataset, fname=None, order=2):
+    ''' Plot the error and its theoretical convergence line. '''
+
+    figure = plt.figure(figsize=(12, 7))
+
+    tot = len(dataset)
+    cols = 2
+    rows = tot // 2
+
+    if tot % cols != 0:
+        rows += 1
+
+    for k, data in enumerate(dataset):
+
+        ax = figure.add_subplot(rows, cols, k+1)
+        ax.set(
+            title=f"{data[0]}",
+            xlabel="$N$",
+            ylabel="L2 Norm of Error",
+            yscale="log",
+            xscale="log"
+        )
+
+        err = np.array(data[1])
+        nx = np.array(data[2])
+
+        ax.scatter(nx, err, marker='x', color='k', label="Solutions")
+        ax.plot(nx, err[-1]*(nx[-1]/nx)**2, linestyle="--",
+                label=r"$\mathcal{O}$" + fr"$(\Delta x^{order})$")
+
+        ax.legend()
+
+    plt.tight_layout()
+
+    if fname is not None:
+        plt.savefig(fname, format="pdf", bbox_inches="tight")
+
+    plt.show()
+
+
+def main():
+    ''' main driver '''
+
+    if len(args.input_file) < 3:
+        raise ValueError('''Must use at least 3 plotfiles with 3 different
+        resolutions that differ by the same multiplicative factor''')
+
+    field_names = ["Variable Name", "L2 Norm(1) (Finest)"]
+
+    for i in range(len(args.input_file[2:])):
+        field_names.append(f"Order of Conv({i+1}:{i+2})")
+        field_names.append(f"L2 Norm({i+2})")
+
+    field_names[-1] += " (Coarsest)"
+
+    table = PrettyTable(field_names)
+
+    data_file = []
+
+    temp_file = io.read(args.input_file[0])
+    for variable in temp_file.cc_data.names:
+
+        l2norms = []
+        nx_list = []
+        row = [variable]
+
+        fdata = io.read(args.input_file[0])
+        cdata = io.read(args.input_file[1])
+        _, l2norm = convergence.compare(fdata.cc_data, cdata.cc_data, variable, args.resolution)
+
+        l2norms.append(l2norm)
+        nx_list.append(cdata.cc_data.grid.nx)
+        row.append(l2norm)
+
+        for n in range(1, len(args.input_file[:-1])):
+            fdata = io.read(args.input_file[n])
+            cdata = io.read(args.input_file[n+1])
+
+            _, l2norm = convergence.compare(fdata.cc_data, cdata.cc_data, variable, args.resolution)
+
+            order_conv = np.sqrt(l2norm/l2norms[-1])
+
+            l2norms.append(l2norm)
+            nx_list.append(cdata.cc_data.grid.nx)
+            row.extend([order_conv, l2norm])
+
+        data_file.append([variable, l2norms, nx_list])
+        table.add_row(row)
+
+    print(table)
+    convergence_plot(data_file, fname=args.out, order=args.order)
+
+
+if __name__ == "__main__":
+    main()

requirements.txt

@@ -15,3 +15,4 @@ sphinxcontrib_bibtex
 sphinx-copybutton
 sphinx-prompt
 importlib-metadata
+prettytable