# Step P1: Import needed NRPy+ core modules:
from outputC import lhrh,outCfunction # NRPy+: Core C code output module
import finite_difference as fin       # NRPy+: Finite difference C code generation module
import NRPy_param_funcs as par        # NRPy+: Parameter interface
import grid as gri                    # NRPy+: Functions having to do with numerical grids
import reference_metric as rfm        # NRPy+: Reference metric support
import deprecated_reference_metric as evil_rfm  # NRPy+: PLEASE DON'T USE.
import cmdline_helper as cmd          # NRPy+: Multi-platform Python command-line interface
import shutil, os, sys                # Standard Python modules for multiplatform OS-level functions

# Step P2: Create C code output directory:
Ccodesdir = os.path.join("ScalarWaveCurvilinear_Playground_Ccodes/")
# First remove C code output directory if it exists
# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty
shutil.rmtree(Ccodesdir, ignore_errors=True)
# Then create a fresh directory
cmd.mkdir(Ccodesdir)

# Step P3: Create executable output directory:
outdir = os.path.join(Ccodesdir,"output/")
cmd.mkdir(outdir)

# Step 1: Set the spatial dimension parameter
#         to three this time, and then read
#         the parameter as DIM.
par.set_parval_from_str("grid::DIM",3)
DIM = par.parval_from_str("grid::DIM")

# Step 2: Set some core parameters, including CoordSystem, boundary condition,
#                                             MoL, timestepping algorithm, FD order,
#                                             floating point precision, and CFL factor:

# Step 2.a: Set the coordinate system for the numerical grid
# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,
#              SymTP, SinhSymTP
CoordSystem     = "SinhSpherical"
par.set_parval_from_str("reference_metric::CoordSystem",CoordSystem)
rfm.reference_metric()

# Step 2.b: Set defaults for Coordinate system parameters.
#           These are perhaps the most commonly adjusted parameters,
#           so we enable modifications at this high level.

# domain_size sets the default value for:
#   * Spherical's params.RMAX
#   * SinhSpherical*'s params.AMAX
#   * Cartesians*'s -params.{x,y,z}min & .{x,y,z}max
#   * Cylindrical's -params.ZMIN & .{Z,RHO}MAX
#   * SinhCylindrical's params.AMPL{RHO,Z}
#   * *SymTP's params.AMAX
domain_size     = 10.0 # Needed for all coordinate systems.

# sinh_width sets the default value for:
#   * SinhSpherical's params.SINHW
#   * SinhCylindrical's params.SINHW{RHO,Z}
#   * SinhSymTP's params.SINHWAA
sinh_width      = 0.4 # If Sinh* coordinates chosen

# sinhv2_const_dr sets the default value for:
#   * SinhSphericalv2's params.const_dr
#   * SinhCylindricalv2's params.const_d{rho,z}
sinhv2_const_dr = 0.05# If Sinh*v2 coordinates chosen

# SymTP_bScale sets the default value for:
#   * SinhSymTP's params.bScale
SymTP_bScale    = 1.0 # If SymTP chosen

# Step 2.c: Set the order of spatial and temporal derivatives;
#           the core data type, and the CFL factor.
# RK_method choices include: Euler, "RK2 Heun", "RK2 MP", "RK2 Ralston", RK3, "RK3 Heun", "RK3 Ralston",
#              SSPRK3, RK4, DP5, DP5alt, CK5, DP6, L6, DP8
RK_method = "RK4"
FD_order  = 4        # Finite difference order: even numbers only, starting with 2. 12 is generally unstable
REAL      = "double" # Best to use double here.
CFL_FACTOR= 1.0

# Step 3: Generate Runge-Kutta-based (RK-based) timestepping code.
#         Each RK substep involves two function calls:
#      3.A: Evaluate RHSs (RHS_string)
#      3.B: Apply boundary conditions (post_RHS_string)
import MoLtimestepping.C_Code_Generation as MoL
from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict
RK_order  = Butcher_dict[RK_method][1]
cmd.mkdir(os.path.join(Ccodesdir,"MoLtimestepping/"))

RHS_string = "rhs_eval(&rfmstruct, &params, RK_INPUT_GFS, RK_OUTPUT_GFS);"

post_RHS_string = "apply_bcs_curvilinear(&params, &bcstruct, NUM_EVOL_GFS, evol_gf_parity, RK_OUTPUT_GFS);"

MoL.MoL_C_Code_Generation(RK_method, RHS_string = RHS_string, post_RHS_string = post_RHS_string,
                          outdir = os.path.join(Ccodesdir,"MoLtimestepping/"))

# Step 4: Import the ScalarWave.InitialData module.
#         This command only declares ScalarWave initial data
#         parameters and the InitialData() function.
import ScalarWave.InitialData as swid

# Step 5: Import ScalarWave_RHSs module.
#         This command only declares ScalarWave RHS parameters
#         and the ScalarWave_RHSs function (called later)
import ScalarWave.ScalarWaveCurvilinear_RHSs as swrhs

# Step 6: Set the finite differencing order to FD_order (set above).
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER",FD_order)

# Step 7: Call the InitialData() function to set up initial data.
#         Options include:
#    "PlaneWave": monochromatic (single frequency/wavelength) plane wave
#    "SphericalGaussian": spherically symmetric Gaussian, with default stdev=3
swid.InitialData(CoordSystem=CoordSystem,Type="PlaneWave")

# Step 8: Generate SymPy symbolic expressions for
#         uu_rhs and vv_rhs; the ScalarWave RHSs.
#         This function also declares the uu and vv
#         gridfunctions, which need to be declared
#         to output even the initial data to C file.
cmd.mkdir(os.path.join(Ccodesdir,"rfm_files/"))
par.set_parval_from_str("reference_metric::enable_rfm_precompute","True")
par.set_parval_from_str("reference_metric::rfm_precompute_Ccode_outdir",os.path.join(Ccodesdir,"rfm_files/"))
swrhs.ScalarWaveCurvilinear_RHSs()
# Step 8.a: Now that we are finished with all the rfm hatted
#           quantities, let's restore them to their closed-
#           form expressions.
par.set_parval_from_str("reference_metric::enable_rfm_precompute","False") # Reset to False to disable rfm_precompute.
rfm.ref_metric__hatted_quantities()

# Step 9: Copy SIMD/SIMD_intrinsics.h to $Ccodesdir/SIMD/SIMD_intrinsics.h
cmd.mkdir(os.path.join(Ccodesdir,"SIMD"))
shutil.copy(os.path.join("SIMD/")+"SIMD_intrinsics.h",os.path.join(Ccodesdir,"SIMD/"))

# Step 10: Generate all needed C functions
enable_FD_functions = False
par.set_parval_from_str("finite_difference::enable_FD_functions",enable_FD_functions)

desc="Part P3: Declare the function for the exact solution at a single point. time==0 corresponds to the initial data."
name="exact_solution_single_point"
outCfunction(
    outfile  = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
    params   ="const REAL xx0,const REAL xx1,const REAL xx2,const paramstruct *restrict params,REAL *uu_exact,REAL *vv_exact",
    body     = fin.FD_outputC("returnstring",[lhrh(lhs="*uu_exact",rhs=swid.uu_ID),
                                              lhrh(lhs="*vv_exact",rhs=swid.vv_ID)]),
    loopopts = "")

desc="Part P4: Declare the function for the exact solution at all points. time==0 corresponds to the initial data."
name="exact_solution_all_points"
outCfunction(
    outfile  = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *restrict params,REAL *restrict xx[3], REAL *restrict in_gfs",
    body     ="""exact_solution_single_point(xx[0][i0],xx[1][i1],xx[2][i2],params,
                            &in_gfs[IDX4S(UUGF,i0,i1,i2)],&in_gfs[IDX4S(VVGF,i0,i1,i2)]);""",
    loopopts = "AllPoints")

desc="Part P5: Declare the function to evaluate the scalar wave RHSs"
includes = None
if enable_FD_functions:
    includes = ["finite_difference_functions.h"]
name="rhs_eval"
outCfunction(
    outfile  = os.path.join(Ccodesdir,name+".h"), includes=includes, desc=desc, name=name,
    params   ="""rfm_struct *restrict rfmstruct,const paramstruct *restrict params,
                 const REAL *restrict in_gfs, REAL *restrict rhs_gfs""",
    body     =fin.FD_outputC("returnstring",[lhrh(lhs=gri.gfaccess("rhs_gfs","uu"),rhs=swrhs.uu_rhs),
                                             lhrh(lhs=gri.gfaccess("rhs_gfs","vv"),rhs=swrhs.vv_rhs)],
                             params="enable_SIMD=True"),
    loopopts = "InteriorPoints,enable_SIMD,enable_rfm_precompute")

# Step 10.b Output functions for computing all finite-difference stencils
if enable_FD_functions:
    fin.output_finite_difference_functions_h(path=Ccodesdir)

import CurviBoundaryConditions.CurviBoundaryConditions as cbcs
cbcs.Set_up_CurviBoundaryConditions(os.path.join(Ccodesdir,"boundary_conditions/"),
                                    Cparamspath=os.path.join("../"))

%%writefile $Ccodesdir/ScalarWaveCurvilinear_Playground.c

// Step P0: Define REAL and NGHOSTS; and declare CFL_FACTOR. This header is generated in NRPy+.
#include "ScalarWaveCurvilinear_Playground_REAL__NGHOSTS__CFL_FACTOR.h"

#include "rfm_files/rfm_struct__declare.h"

#include "declare_Cparameters_struct.h"

// All SIMD intrinsics used in SIMD-enabled C code loops are defined here:
#include "SIMD/SIMD_intrinsics.h"

// Step P1: Import needed header files
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "stdint.h" // Needed for Windows GCC 6.x compatibility
#ifndef M_PI
#define M_PI 3.141592653589793238462643383279502884L
#endif
#ifndef M_SQRT1_2
#define M_SQRT1_2 0.707106781186547524400844362104849039L
#endif

// Step P2: Declare the IDX4S(gf,i,j,k) macro, which enables us to store 4-dimensions of
//           data in a 1D array. In this case, consecutive values of "i"
//           (all other indices held to a fixed value) are consecutive in memory, where
//           consecutive values of "j" (fixing all other indices) are separated by
//           Nxx_plus_2NGHOSTS0 elements in memory. Similarly, consecutive values of
//           "k" are separated by Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1 in memory, etc.
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )
#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) ) ) )
#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \
  for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++)
#define LOOP_ALL_GFS_GPS(ii) _Pragma("omp parallel for") \
  for(int (ii)=0;(ii)<Nxx_plus_2NGHOSTS_tot*NUM_EVOL_GFS;(ii)++)

// Step P3: Set UUGF and VVGF macros, as well as xx_to_Cart()
#include "boundary_conditions/gridfunction_defines.h"

// Step P4: Set xx_to_Cart(const paramstruct *restrict params,
//                     REAL *restrict xx[3],
//                     const int i0,const int i1,const int i2,
//                     REAL xCart[3]),
//           which maps xx->Cartesian via
//    {xx[0][i0],xx[1][i1],xx[2][i2]}->{xCart[0],xCart[1],xCart[2]}
#include "xx_to_Cart.h"

// Step P5: Defines set_Nxx_dxx_invdx_params__and__xx(const int EigenCoord, const int Nxx[3],
//                                       paramstruct *restrict params, REAL *restrict xx[3]),
//          which sets params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for
//          the chosen Eigen-CoordSystem if EigenCoord==1, or
//          CoordSystem if EigenCoord==0.
#include "set_Nxx_dxx_invdx_params__and__xx.h"

// Step P6: Include basic functions needed to impose curvilinear
//          parity and boundary conditions.
#include "boundary_conditions/CurviBC_include_Cfunctions.h"

// Step P7: Find the CFL-constrained timestep
#include "find_timestep.h"

// Part P8: Declare the function for the exact solution at a single point. time==0 corresponds to the initial data.
#include "exact_solution_single_point.h"

// Part P9: Declare the function for the exact solution at all points. time==0 corresponds to the initial data.
#include "exact_solution_all_points.h"

// Part P10: Declare the function to evaluate the scalar wave RHSs
#include "rhs_eval.h"

// main() function:
// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates
// Step 1: Set up scalar wave initial data
// Step 2: Output relative error between numerical and exact solution.
// Step 3: Evolve scalar wave initial data forward in time using Method of Lines with chosen RK-like algorithm,
//         applying quadratic extrapolation outer boundary conditions.
// Step 4: Free all allocated memory
int main(int argc, const char *argv[]) {
    paramstruct params;
#include "set_Cparameters_default.h"

    // Step 0a: Read command-line input, error out if nonconformant
    if(argc != 4 || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < NGHOSTS || atoi(argv[3]) < NGHOSTS) {
        printf("Error: Expected one command-line argument: ./ScalarWaveCurvilinear_Playground Nx0 Nx1 Nx2,\n");
        printf("where Nx[0,1,2] is the number of grid points in the 0, 1, and 2 directions.\n");
        printf("Nx[] MUST BE larger than NGHOSTS (= %d)\n",NGHOSTS);
        exit(1);
    }
    // Step 0b: Set up numerical grid structure, first in space...
    const int Nxx[3] = { atoi(argv[1]), atoi(argv[2]), atoi(argv[3]) };
    if(Nxx[0]%2 != 0 || Nxx[1]%2 != 0 || Nxx[2]%2 != 0) {
        printf("Error: Cannot guarantee a proper cell-centered grid if number of grid cells not set to even number.\n");
        printf("       For example, in case of angular directions, proper symmetry zones will not exist.\n");
        exit(1);
    }

    // Step 0c: Set free parameters, overwriting Cparameters defaults
    //          by hand or with command-line input, as desired.
#include "free_parameters.h"

   // Step 0d: Uniform coordinate grids are stored to *xx[3]
    REAL *xx[3];
    // Step 0d.i: Set bcstruct
    bc_struct bcstruct;
    {
        int EigenCoord = 1;
        // Step 0d.ii: Call set_Nxx_dxx_invdx_params__and__xx(), which sets
        //             params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the
        //             chosen Eigen-CoordSystem.
        set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &params, xx);
        // Step 0d.iii: Set Nxx_plus_2NGHOSTS_tot
#include "set_Cparameters-nopointer.h"
        const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;
        // Step 0e: Find ghostzone mappings; set up bcstruct
#include "boundary_conditions/driver_bcstruct.h"
        // Step 0e.i: Free allocated space for xx[][] array
        for(int i=0;i<3;i++) free(xx[i]);
    }

    // Step 0f: Call set_Nxx_dxx_invdx_params__and__xx(), which sets
    //          params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the
    //          chosen (non-Eigen) CoordSystem.
    int EigenCoord = 0;
    set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &params, xx);

    // Step 0g: Set all C parameters "blah" for params.blah, including
    //          Nxx_plus_2NGHOSTS0 = params.Nxx_plus_2NGHOSTS0, etc.
#include "set_Cparameters-nopointer.h"
    const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;

    // Step 0h: Time coordinate parameters
    const REAL t_final =  0.7*domain_size; /* Final time is set so that at t=t_final,
                                            * data at the origin have not been corrupted
                                            * by the approximate outer boundary condition */

    // Step 0i: Set timestep based on smallest proper distance between gridpoints and CFL factor
    REAL dt = find_timestep(&params, xx);
    //printf("# Timestep set to = %e\n",(double)dt);
    int N_final = (int)(t_final / dt + 0.5); // The number of points in time.
                                             // Add 0.5 to account for C rounding down
                                             // typecasts to integers.
    int output_every_N = (int)((REAL)N_final/800.0);
    if(output_every_N == 0) output_every_N = 1;

    // Step 0j: Error out if the number of auxiliary gridfunctions outnumber evolved gridfunctions.
    //              This is a limitation of the RK method. You are always welcome to declare & allocate
    //              additional gridfunctions by hand.
    if(NUM_AUX_GFS > NUM_EVOL_GFS) {
        printf("Error: NUM_AUX_GFS > NUM_EVOL_GFS. Either reduce the number of auxiliary gridfunctions,\n");
        printf("       or allocate (malloc) by hand storage for *diagnostic_output_gfs. \n");
        exit(1);
    }

    // Step 0k: Allocate memory for gridfunctions
#include "MoLtimestepping/RK_Allocate_Memory.h"

    // Step 0l: Set up precomputed reference metric arrays
    // Step 0l.i: Allocate space for precomputed reference metric arrays.
#include "rfm_files/rfm_struct__malloc.h"

    // Step 0l.ii: Define precomputed reference metric arrays.
    {
#include "set_Cparameters-nopointer.h"
#include "rfm_files/rfm_struct__define.h"
    }

    // Step 1: Set up initial data to be exact solution at time=0:
    params.time = 0.0; exact_solution_all_points(&params, xx, y_n_gfs);

    for(int n=0;n<=N_final;n++)
       { // Main loop to progress forward in time.

        // Step 1a: Set current time to correct value & compute exact solution
        params.time = ((REAL)n)*dt;

        // Step 2: Code validation: Compute log of L2 norm of difference
        //         between numerical and exact solutions:
        //   log_L2_Norm = log10( sqrt[Integral( [numerical - exact]^2 * dV)] ),
        //         where integral is within 30% of the grid outer boundary (domain_size)
        if(n%output_every_N == 0) {
            REAL integral = 0.0;
            REAL numpts   = 0.0;
#pragma omp parallel for reduction(+:integral,numpts)
            LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,
                        NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,
                        NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS) {
                REAL xCart[3]; xx_to_Cart(&params,xx,i0,i1,i2, xCart);
                if(sqrt(xCart[0]*xCart[0] + xCart[1]*xCart[1] + xCart[2]*xCart[2]) < domain_size*0.3) {
                    REAL uu_exact,vv_exact; exact_solution_single_point(xx[0][i0],xx[1][i1],xx[2][i2],&params,
                                                                        &uu_exact,&vv_exact);
                    double num   = (double)y_n_gfs[IDX4S(UUGF,i0,i1,i2)];
                    double exact = (double)uu_exact;
                    integral += (num - exact)*(num - exact);
                    numpts   += 1.0;
                }
            }
            // Compute and output the log of the L2 norm.
            REAL log_L2_Norm = log10(sqrt(integral/numpts));
            printf("%e %e\n",(double)params.time,log_L2_Norm);
        }

        // Step 3: Step forward one timestep (t -> t+dt) in time using
        //           chosen RK-like MoL timestepping algorithm
#include "MoLtimestepping/RK_MoL.h"

    } // End main loop to progress forward in time.

    // Step 4: Free all allocated memory
#include "rfm_files/rfm_struct__freemem.h"
#include "boundary_conditions/bcstruct_freemem.h"
#include "MoLtimestepping/RK_Free_Memory.h"
    for(int i=0;i<3;i++) free(xx[i]);
    return 0;
}

import cmdline_helper as cmd

cmd.C_compile(os.path.join(Ccodesdir,"ScalarWaveCurvilinear_Playground.c"),
              os.path.join(outdir,"ScalarWaveCurvilinear_Playground"),compile_mode="optimized")
# !clang -Ofast -fopenmp -mavx2 -mfma ScalarWave/ScalarWaveCurvilinear_Playground.c -o ScalarWaveCurvilinear_Playground -lm
# !icc -align -qopenmp -xHost -O2 -qopt-report=5 -qopt-report-phase ipo -qopt-report-phase vec -vec-threshold1 -qopt-prefetch=4 ScalarWave/ScalarWaveCurvilinear_Playground.c -o ScalarWaveCurvilinear_Playground
# !gcc-7 -Ofast -fopenmp -march=native ScalarWave/ScalarWaveCurvilinear_Playground.c -o ScalarWaveCurvilinear_Playground -lm

# Change to output directory
os.chdir(outdir)
# Clean up existing output files
cmd.delete_existing_files("out-*resolution.txt")
# Run executable
if par.parval_from_str("reference_metric::CoordSystem") == "Cartesian":
    cmd.Execute("ScalarWaveCurvilinear_Playground", "16 16 16", "out-lowresolution.txt")
    cmd.Execute("ScalarWaveCurvilinear_Playground", "24 24 24", "out-medresolution.txt")
else:
    cmd.Execute("ScalarWaveCurvilinear_Playground", "16  8 16", "out-lowresolution.txt")
    # 4.28s with icc and FD order = 10.
    cmd.Execute("ScalarWaveCurvilinear_Playground", "24 12 24", "out-medresolution.txt")
    ########################################
    # BENCHMARK 48x24x48 RUN, FD order = 4. desktop: 17.33s
    # laptop: 51.82s on icc. 45.02s on GCC 9, 45.03s on GCC 7, 51.67s on clang
    #     cmd.Execute("ScalarWaveCurvilinear_Playground", "48 24 48", "out-hghresolution.txt")

# %timeit cmd.Execute("ScalarWaveCurvilinear_Playground", "48 24 48", "out-hghresolution.txt", verbose=False)
# FD functions disabled:
# 16 s ± 702 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
# FD functions enabled:
# 16.1 s ± 384 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

# Return to root directory
os.chdir(os.path.join("../../"))

%matplotlib inline
import matplotlib.pyplot as plt
import mpmath as mp
import csv

def file_reader(filename):
    with open(filename) as file:
        reader = csv.reader(file, delimiter=" ")
        data  = list(zip(*reader))
        # data is a tuple of strings. Tuples are immutable, and we need to perform math on
        #   the data, so here we convert tuple to lists of floats:
        data0 = []
        data1 = []
        for i in range(len(data[0])):
            data0.append(float(data[0][i]))
            data1.append(float(data[1][i]))
        return data0,data1

first_col16,second_col16 = file_reader(os.path.join(outdir,'out-lowresolution.txt'))
first_col24,second_col24 = file_reader(os.path.join(outdir,'out-medresolution.txt'))

second_col16_rescaled4o = []
second_col16_rescaled5o = []
for i in range(len(second_col16)):
    # data16 = data24*(16/24)**4
    # -> log10(data24) = log10(data24) + 4*log10(16/24)
    second_col16_rescaled4o.append(second_col16[i] + 4*mp.log10(16./24.))
    second_col16_rescaled5o.append(second_col16[i] + 5*mp.log10(16./24.))

# https://matplotlib.org/gallery/text_labels_and_annotations/legend.html#sphx-glr-gallery-text-labels-and-annotations-legend-py
fig, ax = plt.subplots()

plt.title("Demonstrating 4th-order Convergence: "+par.parval_from_str("reference_metric::CoordSystem")+" Coordinates")
plt.xlabel("time")
plt.ylabel("log10(Max relative error)")

ax.plot(first_col24, second_col24, 'k-', label='logErel(N0=24)')
ax.plot(first_col16, second_col16_rescaled4o, 'k--', label='logErel(N0=16) + log((16/24)^4)')
ax.set_ylim([-8.05,-1.7]) # Manually set the y-axis range case, since the log10
                          # relative error at t=0 could be -inf or about -16,
                          # resulting in very different-looking plots
                          # despite the data being the same to roundoff.
if par.parval_from_str("reference_metric::CoordSystem") == "Cartesian":
    ax.set_ylim([-2.68,-1.62])
if par.parval_from_str("reference_metric::CoordSystem") == "Cylindrical":
    ax.plot(first_col16, second_col16_rescaled5o, 'k.',  label='(Assuming 5th-order convergence)')
legend = ax.legend(loc='lower right', shadow=True, fontsize='large')
legend.get_frame().set_facecolor('C1')
plt.show()

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Start_to_Finish-ScalarWaveCurvilinear")