#!/usr/bin/env python
# coding: utf-8

# <script async src="https://www.googletagmanager.com/gtag/js?id=UA-59152712-8"></script>
# <script>
#   window.dataLayer = window.dataLayer || [];
#   function gtag(){dataLayer.push(arguments);}
#   gtag('js', new Date());
# 
#   gtag('config', 'UA-59152712-8');
# </script>
# 
# # `MaxwellVacuum`: Solving Maxwell Equations in Vacuum with NRPy+ and the Einstein Toolkit
# 
# ## Authors: Terrence Pierre Jacques, Patrick Nelson, & Zach Etienne
# 
# ## This notebook generates `MaxwellVacuum`, an [Einstein Toolkit](https://einsteintoolkit.org) thorn for solving Maxwell's equations in vacuum, in flat spacetime and Cartesian coordinates. This thorn is highly optimized for modern CPU architectures, featuring SIMD intrinsics and OpenMP support.
# 
# ### Formatting improvements courtesy Brandon Clark
# 
# [comment]: <> (Abstract: TODO)
# 
# **Notebook Status:** <font color='green'><b> Validated</b></font> 
# 
# **Validation Notes:** As demonstrated [in the plot below](#code_validation), the numerical errors converge to zero at the expected rate, and we observe similar qualitatitve behavior of of the error nodes for Systems I (ADM - like) and II (BSSN - like). Also validated against the NRPy+ [Tutorial-Start_to_Finish-Solving_Maxwells_Equations_in_Vacuum-Cartesian](Tutorial-Start_to_Finish-Solving_Maxwells_Equations_in_Vacuum-Cartesian.ipynb) notebook.
# 
# ### NRPy+ Source Code for this module: 
# * [Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py](../edit/Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py) [\[**tutorial**\]](Tutorial-VacuumMaxwell_formulation_Cartesian.ipynb)
# 
# ## Introduction:
# This tutorial notebook constructs an Einstein Toolkit (ETK) thorn (module) that will set up expressions for the right-hand sides of Maxwell's equations as described in [Knapp, Walker, & Baumgarte (2002)](https://arxiv.org/abs/gr-qc/0201051) and the [NRPy+ `Tutorial-VacuumMaxwell_formulation_Cartesian` Jupyter notebook](Tutorial-VacuumMaxwell_formulation_Cartesian.ipynb), both for System I:
# 
# \begin{align}
# \partial_t A^i &= -E^i - \partial^i \varphi, \\
# \partial_t E^i &= \partial^i \partial_j A^j - \partial_j \partial^j A^i, \\
# \partial_t \varphi &= -\partial_i A^i,
# \end{align}
# 
# with associated constraint
# 
# $$
# \partial_i E^i = 0,
# $$
# 
# and for System II:
# 
# \begin{align}
# \partial_t A^i &= -E^i - \partial^i \varphi, \\
# \partial_t E^i &= \partial^i \Gamma - \partial_j \partial^j A^i, \\
# \partial_t \Gamma &= - \partial_i \partial^i \varphi, \\
# \partial_t \varphi &= -\Gamma,
# \end{align}
# 
# with associated constraints
# 
# \begin{align}
# \Gamma - \partial_i A^i &= 0, \\
# \partial_i E^i &= 0.
# \end{align}
# 
# This thorn is largely based on and should function similarly to the [`BaikalETK`](Tutorial-BaikalETK.ipynb) thorns, which solve Einstein's equations of general relativity in Cartesian coordinates, in the BSSN formalism. Further, we generate the C kernels for a number of finite difference orders, so that users are free to choose the finite-differencing order at runtime.
# 
# When interfaced properly with the ETK, this module will propagate the initial data for $E^i$, $A^i$, and $\varphi$ defined in [`MaxwellVacuumID`](Tutorial-ETK_thorn-MaxwellVacuumID.ipynb) forward in time by integrating the equations for $\partial_t E^i$, $\partial_t A^i$ and $\partial_t \varphi$ subject to spatial boundary conditions. The time evolution itself is handled by the `MoL` (Method of Lines) thorn in the `CactusNumerical` arrangement, and the boundary conditions by the `Boundary` thorn in the `CactusBase` arrangement. Specifically, we implement ETK's `NewRad` boundary condition driver, i.e. a radiation (Sommerfeld) boundary condition. 
# 
# Similar to the [`MaxwellVacuumID`](Tutorial-ETK_thorn-MaxwellVacuumID.ipynb) module, we will construct the MaxwellVacuum module in two steps.
# 
# 1. Call on NRPy+ to convert the SymPy expressions for the evolution equations into C-code kernels.
# 1. Write the C code drivers and linkages to the Einstein Toolkit infrastructure (i.e., the .ccl files) to complete this Einstein Toolkit module.
# 
# The above steps will be performed in this notebook by first outputting *Python* code into a module, and then loading the Python module to generate the C code *in parallel* (using the [`multiprocessing`](https://docs.python.org/3/library/multiprocessing.html) built-in Python module).

# <a id='toc'></a>
# 
# # Table of Contents
# $$\label{toc}$$
# 
# This notebook is organized as follows
# 
# 1. [Step 1](#initializenrpy): Initialize needed Python/NRPy+ modules
# 1. [Step 2](#mwev): NRPy+-generated C code kernels for Maxwell spacetime solve
#     1. [Step 2.a](#mwevrhs): Maxwell RHS expressions
#     1. [Step 2.b](#constraints): Divergence Constraint
#     1. [Step 2.c](#driver_call_codegen_funcs):  Given `WhichPart` parameter choice, generate algorithm for calling corresponding function within `MaxwellVacuum_C_kernels_codegen_onepart()` to generate C code kernel 
#     1. [Step 2.e](#kernel_codegen): Generate C code kernels for `MaxwellVacuum`
#         1. [Step 2.e.i](#feature_choice): Set compile-time and runtime parameters for `MaxwellVacuum`
#         1. [Step 2.e.ii](#parallel_codegen): Generate all C-code kernels for `MaxwellVacuum`, in parallel if possible
# 1. [Step 3](#cclfiles): CCL files - Define how this module interacts and interfaces with the wider Einstein Toolkit infrastructure
#     1. [Step 3.a](#paramccl): `param.ccl`: specify free parameters within `MaxwellVacuum`
#     1. [Step 3.b](#interfaceccl): `interface.ccl`: define needed gridfunctions; provide keywords denoting what this thorn provides and what it should inherit from other thorns
#     1. [Step 3.c](#scheduleccl): `schedule.ccl`:schedule all functions used within `MaxwellVacuum`, specify data dependencies within said functions, and allocate memory for gridfunctions
# 1. [Step 4](#cdrivers): C driver functions for ETK registration & NRPy+ generated kernels
#     1. [Step 4.a](#etkfunctions): Needed ETK functions: Banner, Symmetry registration, Parameter sanity check,  Method of Lines (`MoL`) registration, Boundary condition
#     1. [Step 4.b](#mwevrhss) Evaluate Maxwell right-hand-sides (RHSs)
#     1. [Step 4.c](#diagnostics): Diagnostics: Computing the divergence constraint
# 
#     1. [Step 4.d](#outcdrivers): Output all C driver functions needed for `MaxwellVacuum`
#     1. [Step 4.e](#makecodedefn): `make.code.defn`: List of all C driver functions needed to compile `MaxwellVacuum`
# 1. [Step 5](#code_validation): Code Validation
#     1. [Step 5.a](#convergence): Error Convergence
#     1. [Step 5.b](#errornodes): Behavior of Error Nodes
# 1. [Step 6](#latex_pdf_output): Output this notebook to $\LaTeX$-formatted PDF file

# <a id='initializenrpy'></a>
# 
# # Step 1: Initialize needed Python/NRPy+ modules \[Back to [top](#toc)\]
# 
# $$\label{initializenrpy}$$

# In[1]:


MaxwellVacuumdir = "MaxwellNRPy_py_dir"
import cmdline_helper as cmd  # NRPy+: Multi-platform Python command-line interface
import os, sys                # Standard Python modules for multiplatform OS-level functions

cmd.mkdir(os.path.join(MaxwellVacuumdir))
# Write an empty __init__.py file in this directory so that Python2 can load modules from it.
with open(os.path.join(MaxwellVacuumdir, "__init__.py"), "w") as file:
    pass


# In[2]:


get_ipython().run_cell_magic('writefile', '$MaxwellVacuumdir/MaxwellVacuum_C_kernels_codegen.py', '\n# Step 1: Import needed core NRPy+ modules\nfrom outputC import lhrh         # NRPy+: Core C code output module\nimport finite_difference as fin  # NRPy+: Finite difference C code generation module\nimport NRPy_param_funcs as par   # NRPy+: Parameter interface\nimport grid as gri               # NRPy+: Functions having to do with numerical grids\nimport loop as lp                # NRPy+: Generate C code loops\nimport indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\nimport reference_metric as rfm   # NRPy+: Reference metric support\nimport os, sys                   # Standard Python modules for multiplatform OS-level functions\nimport time                      # Standard Python module; useful for benchmarking\nimport Maxwell.VacuumMaxwell_Flat_Evol_Cartesian as rhs\n\ndef MaxwellVacuum_C_kernels_codegen_onepart(params=\n        "WhichPart=Maxwell_RHSs,ThornName=MaxwellVacuum,FD_order=4"):\n    # Set default parameters\n    WhichPart = "Maxwell_RHSs"\n    ThornName = "MaxwellVacuum"\n    FD_order = 4\n\n    import re\n    if params != "":\n        params2 = re.sub("^,","",params)\n        params = params2.strip()\n        splitstring = re.split("=|,", params)\n\n#         if len(splitstring) % 2 != 0:\n#             print("outputC: Invalid params string: "+params)\n#             sys.exit(1)\n\n        parnm = []\n        value = []\n        for i in range(int(len(splitstring)/2)):\n            parnm.append(splitstring[2*i])\n            value.append(splitstring[2*i+1])\n\n        for i in range(len(parnm)):\n            parnm.append(splitstring[2*i])\n            value.append(splitstring[2*i+1])\n\n        for i in range(len(parnm)):\n            if parnm[i] == "WhichPart":\n                WhichPart = value[i]\n            elif parnm[i] == "ThornName":\n                ThornName = value[i]\n            elif parnm[i] == "FD_order":\n                FD_order = int(value[i])\n            else:\n                print("MaxwellVacuum Error: Could not parse input param: "+parnm[i])\n                sys.exit(1)\n\n    # Set output directory for C kernels\n    outdir = os.path.join(ThornName, "src")  # Main C code output directory\n\n    # Set spatial dimension (must be 3 for Maxwell)\n    par.set_parval_from_str("grid::DIM",3)\n\n    # Step 2: Set some core parameters, including CoordSystem MoL timestepping algorithm,\n    #                                 FD order, floating point precision, and CFL factor:\n    # Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,\n    #              SymTP, SinhSymTP\n    # NOTE: Only CoordSystem == Cartesian makes sense here; new\n    #       boundary conditions are needed within the ETK for\n    #       Spherical, etc. coordinates.\n    CoordSystem     = "Cartesian"\n\n    par.set_parval_from_str("reference_metric::CoordSystem",CoordSystem)\n    rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating Maxwell RHSs, etc.\n\n    # Set the gridfunction memory access type to ETK-like, so that finite_difference\n    #    knows how to read and write gridfunctions from/to memory.\n    par.set_parval_from_str("grid::GridFuncMemAccess","ETK")\n')


# <a id='mwev'></a>
# 
# # Step 2: NRPy+-generated C code kernels for solving \[Back to [top](#toc)\]
# $$\label{mwev}$$
# 
# <a id='mwevrhs'></a>
# 
# ## Step 2.a: Maxwell RHS expressions \[Back to [top](#toc)\]
# $$\label{mwevrhs}$$
# 
# `MaxwellVacuum` implements Maxwell's equations in flat spacetime and in Cartesian coordinates, which is fully documented within NRPy+ ([start here](Tutorial-VacuumMaxwell_formulation_Cartesian.ipynb)). 
# 
# Here, we simply call the [Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py](../edit/Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py); [\[**tutorial**\]](Tutorial-VacuumMaxwell_formulation_Cartesian.ipynb) NRPy+ Python module to generate the symbolic expressions and then output the finite-difference C code form of the equations using NRPy+'s [finite_difference](finite_difference.py) ([**tutorial**](Tutorial-Finite_Difference_Derivatives.ipynb)) C code generation module.
# 

# In[3]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_kernels_codegen.py', '\n    def Maxwell_RHSs__generate_symbolic_expressions_SystemI():\n        ######################################\n        # START: GENERATE SYMBOLIC EXPRESSIONS\n        print("Generating symbolic expressions for Maxwell RHSs SystemI...\\n")\n        start = time.time()\n\n        # Set which system to use, which are defined in Maxwell/InitialData.py\n        par.initialize_param(par.glb_param("char", "Maxwell.InitialData","System_to_use","System_I"))\n\n        par.set_parval_from_str("reference_metric::enable_rfm_precompute","True")\n        par.set_parval_from_str("reference_metric::rfm_precompute_Ccode_outdir",os.path.join(outdir,"rfm_files/"))\n\n        rhs.VacuumMaxwellRHSs()\n\n        end = time.time()\n        print("(BENCH) Finished Maxwell RHS symbolic expressions in "+str(end-start)+" seconds.\\n")\n        # END: GENERATE SYMBOLIC EXPRESSIONS\n        ######################################\n\n        # Step 2: Register new gridfunctions so they can be written to by NRPy.\n        # System I:\n        AIU = ixp.register_gridfunctions_for_single_rank1("EVOL","AIU")\n        EIU = ixp.register_gridfunctions_for_single_rank1("EVOL","EIU")\n        psiI = gri.register_gridfunctions("EVOL","psiI")\n        C_SystemI = rhs.C\n\n        Maxwell_RHSs_SymbExpressions_SystemI = [\\\n                         lhrh(lhs=gri.gfaccess("rhs_gfs","AIU0"),rhs=rhs.ArhsU[0]),\\\n                         lhrh(lhs=gri.gfaccess("rhs_gfs","AIU1"),rhs=rhs.ArhsU[1]),\\\n                         lhrh(lhs=gri.gfaccess("rhs_gfs","AIU2"),rhs=rhs.ArhsU[2]),\\\n                         lhrh(lhs=gri.gfaccess("rhs_gfs","EIU0"),rhs=rhs.ErhsU[0]),\\\n                         lhrh(lhs=gri.gfaccess("rhs_gfs","EIU1"),rhs=rhs.ErhsU[1]),\\\n                         lhrh(lhs=gri.gfaccess("rhs_gfs","EIU2"),rhs=rhs.ErhsU[2]),\\\n                         lhrh(lhs=gri.gfaccess("rhs_gfs","psi"),rhs=rhs.psi_rhs)]\n\n        # Avoid re-registering GFs later\n        del gri.glb_gridfcs_list[:8]\n\n        return [Maxwell_RHSs_SymbExpressions_SystemI, C_SystemI]\n\n    def Maxwell_RHSs__generate_symbolic_expressions_SystemII():\n        ######################################\n        # START: GENERATE SYMBOLIC EXPRESSIONS\n        print("Generating symbolic expressions for Maxwell RHSs SystemII...\\n")\n        start = time.time()\n\n        # Set which system to use, which are defined in Maxwell/VacuumMaxwell_Flat_Cartesian_ID.py\n        par.set_parval_from_str("Maxwell.InitialData::System_to_use","System_II")\n\n        # Avoid re-registering GFs\n        gri.glb_gridfcs_list = []\n\n        rhs.VacuumMaxwellRHSs()\n\n        end = time.time()\n        print("(BENCH) Finished Maxwell RHS symbolic expressions in "+str(end-start)+" seconds.\\n")\n        # END: GENERATE SYMBOLIC EXPRESSIONS\n        ######################################\n\n        # Register gridfunctions so they can be written to by NRPy.\n        # System II:\n        AIIU = ixp.register_gridfunctions_for_single_rank1("EVOL","AIIU")\n        EIIU = ixp.register_gridfunctions_for_single_rank1("EVOL","EIIU")\n        psiII = gri.register_gridfunctions("EVOL","psiII")\n        GammaII = gri.register_gridfunctions("EVOL","GammaII")\n        C_SystemII = rhs.C\n        G_SystemII = rhs.G\n\n        Maxwell_RHSs_SymbExpressions_SystemII =  [\\\n                          lhrh(lhs=gri.gfaccess("rhs_gfs","AIIU0"),rhs=rhs.ArhsU[0]),\\\n                          lhrh(lhs=gri.gfaccess("rhs_gfs","AIIU1"),rhs=rhs.ArhsU[1]),\\\n                          lhrh(lhs=gri.gfaccess("rhs_gfs","AIIU2"),rhs=rhs.ArhsU[2]),\\\n                          lhrh(lhs=gri.gfaccess("rhs_gfs","EIIU0"),rhs=rhs.ErhsU[0]),\\\n                          lhrh(lhs=gri.gfaccess("rhs_gfs","EIIU1"),rhs=rhs.ErhsU[1]),\\\n                          lhrh(lhs=gri.gfaccess("rhs_gfs","EIIU2"),rhs=rhs.ErhsU[2]),\\\n                          lhrh(lhs=gri.gfaccess("rhs_gfs","psi"),rhs=rhs.psi_rhs),\\\n                          lhrh(lhs=gri.gfaccess("rhs_gfs","Gamma"),rhs=rhs.Gamma_rhs)]\n\n        return [Maxwell_RHSs_SymbExpressions_SystemII, C_SystemII, G_SystemII]\n\n    def Maxwell_RHSs__generate_Ccode(all_RHSs_exprs_list_SystemI, all_RHSs_exprs_list_SystemII):\n        Maxwell_RHSs_SymbExpressions_SystemI  = all_RHSs_exprs_list_SystemI\n        Maxwell_RHSs_SymbExpressions_SystemII = all_RHSs_exprs_list_SystemII\n\n        print("Generating C code for Maxwell RHSs (FD_order="+str(FD_order)+") in "+par.parval_from_str("reference_metric::CoordSystem")+" coordinates.")\n        start = time.time()\n\n        # Store original finite-differencing order:\n        FD_order_orig = par.parval_from_str("finite_difference::FD_CENTDERIVS_ORDER")\n        # Set new finite-differencing order:\n        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order)\n\n        Maxwell_RHSs_string_SystemI = fin.FD_outputC("returnstring",Maxwell_RHSs_SymbExpressions_SystemI,\n                                          params="outCverbose=False,enable_SIMD=True").replace("AU","AIU")\\\n                                                                    .replace("EU","EIU")\\\n                                                                    .replace("psi","psiI")\n\n        Maxwell_RHSs_string_SystemII = fin.FD_outputC("returnstring",Maxwell_RHSs_SymbExpressions_SystemII,\n                                          params="outCverbose=False,enable_SIMD=True").replace("AU","AIIU")\\\n                                                                    .replace("EU","EIIU")\\\n                                                                    .replace("psi","psiII")\\\n                                                                    .replace("Gamma","GammaII")\n\n\n        filename = "Maxwell_RHSs"+"_FD_order_"+str(FD_order)+".h"\n        with open(os.path.join(outdir,filename), "w") as file:\n            file.write(lp.loop(["i2","i1","i0"],["cctk_nghostzones[2]","cctk_nghostzones[1]","cctk_nghostzones[0]"],\n           ["cctk_lsh[2]-cctk_nghostzones[2]","cctk_lsh[1]-cctk_nghostzones[1]","cctk_lsh[0]-cctk_nghostzones[0]"],\n                               ["1","1","SIMD_width"],\n                                ["#pragma omp parallel for",\n                                 "#include \\"rfm_files/rfm_struct__SIMD_outer_read2.h\\"",\n                                 r"""    #include "rfm_files/rfm_struct__SIMD_outer_read1.h"\n    #if (defined __INTEL_COMPILER && __INTEL_COMPILER_BUILD_DATE >= 20180804)\n        #pragma ivdep         // Forces Intel compiler (if Intel compiler used) to ignore certain SIMD vector dependencies\n        #pragma vector always // Forces Intel compiler (if Intel compiler used) to vectorize\n    #endif"""],"",\n                    "#include \\"rfm_files/rfm_struct__SIMD_inner_read0.h\\"\\n"+Maxwell_RHSs_string_SystemI+Maxwell_RHSs_string_SystemII))\n\n        # Restore original finite-differencing order:\n        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order_orig)\n\n        end = time.time()\n        print("(BENCH) Finished Maxwell_RHS C codegen (FD_order="+str(FD_order)+") in " + str(end - start) + " seconds.\\n")\n')


# <a id='constraints'></a>
# 
# ## Step 2.b: Constraints \[Back to [top](#toc)\]
# $$\label{constraints}$$
# 
# Finally output the C code for evaluating the constraints, described in [this tutorial](Tutorial-VacuumMaxwell_Cartesian_RHSs.ipynb). In the absence of numerical error, the constraints should evaluate to zero, but due to numerical (typically truncation and roundoff) error they do not. We will therefore measure constraint violations to gauge the accuracy of our simulation, and ultimately determine whether errors are dominated by numerical finite differencing (truncation) error as expected.

# In[4]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_kernels_codegen.py', '\n    def Maxwell_constraints__generate_symbolic_expressions_and_C_code(C_SystemI,\n                                                                      C_SystemII,\n                                                                      G_SystemII):\n\n        DivEI  = gri.register_gridfunctions("AUX", "DivEI")\n        DivEII = gri.register_gridfunctions("AUX", "DivEII")\n        GII    = gri.register_gridfunctions("AUX", "GII")\n\n        # Store original finite-differencing order:\n        FD_order_orig = par.parval_from_str("finite_difference::FD_CENTDERIVS_ORDER")\n        # Set new finite-differencing order:\n        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order)\n\n        start = time.time()\n        print("Generating optimized C code for constraints.")\n        Constraints_string_SystemI = fin.FD_outputC("returnstring",\n                                        [lhrh(lhs=gri.gfaccess("aux_gfs", "DivEI"),  rhs=C_SystemI)],\n                                        params="outCverbose=False").replace("AU","AIU")\\\n                                                                   .replace("EU","EIU")\\\n                                                                   .replace("psi","psiI")\n\n        Constraints_string_SystemII = fin.FD_outputC("returnstring",\n                                        [lhrh(lhs=gri.gfaccess("aux_gfs", "DivEII"), rhs=C_SystemII),\n                                         lhrh(lhs=gri.gfaccess("aux_gfs", "GII"),    rhs=G_SystemII)],\n                                        params="outCverbose=False").replace("AU","AIIU")\\\n                                                                   .replace("EU","EIIU")\\\n                                                                   .replace("psi","psiII")\\\n                                                                   .replace("Gamma","GammaII")\n\n        with open(os.path.join(outdir,"Maxwell_constraints"+"_FD_order_"+str(FD_order)+".h"), "w") as file:\n            file.write(lp.loop(["i2","i1","i0"],["cctk_nghostzones[2]","cctk_nghostzones[1]","cctk_nghostzones[0]"],\n           ["cctk_lsh[2]-cctk_nghostzones[2]","cctk_lsh[1]-cctk_nghostzones[1]","cctk_lsh[0]-cctk_nghostzones[0]"],\n                               ["1","1","1"],["#pragma omp parallel for","",""], "",Constraints_string_SystemI+Constraints_string_SystemII))\n\n        # Restore original finite-differencing order:\n        par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order_orig)\n\n        end = time.time()\n        print("(BENCH) Finished constraint C codegen (FD_order="+str(FD_order)+") in " + str(end - start) + " seconds.\\n")\n')


# <a id='driver_call_codegen_funcs'></a>
# 
# ## Step 2.d: Given `WhichPart` parameter choice, generate algorithm for calling corresponding function within `MaxwellVacuum_C_kernels_codegen_onepart()` to generate C code kernel \[Back to [top](#toc)\]
# $$\label{driver_call_codegen_funcs}$$

# In[5]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_kernels_codegen.py', '\n\n    SystemI_exprs  = Maxwell_RHSs__generate_symbolic_expressions_SystemI()\n    SystemII_exprs = Maxwell_RHSs__generate_symbolic_expressions_SystemII()\n\n    if WhichPart=="Maxwell_RHSs":\n        Maxwell_RHSs__generate_Ccode(SystemI_exprs[0], SystemII_exprs[0])\n    elif WhichPart=="Maxwell_constraints":\n        Maxwell_constraints__generate_symbolic_expressions_and_C_code(SystemI_exprs[1],\n                                                                     SystemII_exprs[1],\n                                                                     SystemII_exprs[2])\n    else:\n        print("Error: WhichPart = "+WhichPart+" is not recognized.")\n        sys.exit(1)\n\n    # Avoid re-registering GFs\n    del gri.glb_gridfcs_list[7:17]\n\n    # Store all NRPy+ environment variables to an output string so NRPy+ environment from within this subprocess can be easily restored\n    import pickle\n    # https://www.pythonforthelab.com/blog/storing-binary-data-and-serializing/\n    outstr = []\n    outstr.append(pickle.dumps(len(gri.glb_gridfcs_list)))\n    for lst in gri.glb_gridfcs_list:\n        outstr.append(pickle.dumps(lst.gftype))\n        outstr.append(pickle.dumps(lst.name))\n        outstr.append(pickle.dumps(lst.rank))\n        outstr.append(pickle.dumps(lst.DIM))\n\n    outstr.append(pickle.dumps(len(par.glb_params_list)))\n    for lst in par.glb_params_list:\n        outstr.append(pickle.dumps(lst.type))\n        outstr.append(pickle.dumps(lst.module))\n        outstr.append(pickle.dumps(lst.parname))\n        outstr.append(pickle.dumps(lst.defaultval))\n\n    outstr.append(pickle.dumps(len(par.glb_Cparams_list)))\n    for lst in par.glb_Cparams_list:\n        outstr.append(pickle.dumps(lst.type))\n        outstr.append(pickle.dumps(lst.module))\n        outstr.append(pickle.dumps(lst.parname))\n        outstr.append(pickle.dumps(lst.defaultval))\n    return outstr\n')


# <a id='kernel_codegen'></a>
# 
# ## Step 2.c: Generate C code kernels for `MaxwellVacuum` \[Back to [top](#toc)\]
# $$\label{kernel_codegen}$$
# 
# Here we generate the C code kernels (i.e., the C-code representation of the equations needed) for `MaxwellVacuum`.
# 
# <a id='feature_choice'></a>
# 
# ### Step 2.c.i: Set compile-time and runtime parameters for `MaxwellVacuum`  \[Back to [top](#toc)\]
# $$\label{feature_choice}$$
# 
# NRPy+ is a code generation package that is designed to offer maximum flexibility *at the time of C code generation*. As a result, although NRPy+ can in principle output an infinite variety of C code kernels for solving systems of partial differential equations, generally free parameters in each kernel steerable at *runtime* are restricted to simple scalars. This leads to more optimized kernels, but at the expense of flexibility in generating multiple kernels (e.g. one per finite-differencing order). Reducing the number of kernels and adding more flexibility at runtime will be a focus of future work.
# 
# For now, `MaxwellVacuum` supports the following runtime options:
# 
# * `MaxwellVacuum`: Evolution of the Maxwell's equations.
#     * Finite differencing of orders 2, 4, 6, and 8 via runtime parameter `FD_order`
# 
# Next we set up the default parameter lists for `MaxwellVacuum_C_kernels_codegen_onepart()` for the `MaxwellVacuum` thorn. We set these parameter lists as strings to make parallelizing the code generation far easier (easier to pass a list of strings than a list of function arguments to Python's `multiprocessing.Pool()`).

# In[6]:


# Step 2.e.i: Set compile-time and runtime parameters for Maxwell and MaxwellVacuum
#         Runtime parameters for
#         MaxwellVacuum: FD_orders = [2,4,6,8];

paramslist = []
FD_orders = [2,4,6,8]
WhichParamSet = 0
ThornName = "MaxwellVacuum"
for WhichPart in ["Maxwell_RHSs","Maxwell_constraints"]:
    for FD_order in FD_orders:
        paramstr = "WhichPart="+WhichPart+","
        paramstr+= "ThornName="+ThornName+","
        paramstr+= "FD_order="+str(FD_order)+","
        paramslist.append(paramstr)
        WhichParamSet = WhichParamSet + 1

paramslist.sort() # Sort the list alphabetically.


# <a id='parallel_codegen'></a>
# 
# ### Step 2.e.ii: Generate all C-code kernels for `MaxwellVacuum`, in parallel if possible \[Back to [top](#toc)\]
# $$\label{parallel_codegen}$$

# In[7]:


nrpy_dir_path = os.path.join(".")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)

# Create all output directories if they do not yet exist
import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
for ThornName in ["MaxwellVacuum"]:
    outrootdir = ThornName
    cmd.mkdir(os.path.join(outrootdir))
    outdir = os.path.join(outrootdir,"src") # Main C code output directory

    # Copy SIMD/SIMD_intrinsics.h to $outdir/SIMD/SIMD_intrinsics.h, replacing
    #   the line "#define REAL_SIMD_ARRAY REAL" with "#define REAL_SIMD_ARRAY CCTK_REAL"
    #   (since REAL is undefined in the ETK, but CCTK_REAL takes its place)
    cmd.mkdir(os.path.join(outdir,"SIMD"))
    import fileinput
    f = fileinput.input(os.path.join(nrpy_dir_path,"SIMD","SIMD_intrinsics.h"))
    with open(os.path.join(outdir,"SIMD","SIMD_intrinsics.h"),"w") as outfile:
        for line in f:
            outfile.write(line.replace("#define REAL_SIMD_ARRAY REAL", "#define REAL_SIMD_ARRAY CCTK_REAL"))

    # Create directory for rfm_files output
    cmd.mkdir(os.path.join(outdir,"rfm_files"))

# Start parallel C code generation (codegen)
# NRPyEnvVars stores the NRPy+ environment from all the subprocesses in the following
#     parallel codegen
NRPyEnvVars = []

import time   # Standard Python module for benchmarking
import logging
start = time.time()
if __name__ == "__main__":
#     try:
#         if os.name == 'nt':
#             # Windows & Jupyter multiprocessing do not mix, so we run in serial on Windows.
#             #  Here's why: https://stackoverflow.com/questions/45719956/python-multiprocessing-in-jupyter-on-windows-attributeerror-cant-get-attribut
#             raise Exception("Parallel codegen currently not available in Windows")
#         # Step 3.d.ii: Import the multiprocessing module.
#         import multiprocessing
#         print("***************************************")
#         print("Starting parallel C kernel codegen...")
#         print("***************************************")

#         # Step 3.d.iii: Define master function for parallelization.
#         #           Note that lambdifying this doesn't work in Python 3
#         def master_func(i):
#             import MaxwellNRPy_py_dir.MaxwellVacuum_C_kernels_codegen as MCk
#             return MCk.MaxwellVacuum_C_kernels_codegen_onepart(params=paramslist[i])

#         # Step 3.d.iv: Evaluate list of functions in parallel if possible;
#         #           otherwise fallback to serial evaluation:
#         pool = multiprocessing.Pool() #processes=len(paramslist))
#         NRPyEnvVars.append(pool.map(master_func,range(len(paramslist))))
#         pool.terminate()
#         pool.join()
#     except:
#         logging.exception("Ignore this warning/backtrace if on a system in which serial codegen is necessary:")
    print("***************************************")
    print("Starting serial C kernel codegen...")
    print("***************************************")
    # Steps 3.d.ii-iv, alternate: As fallback, evaluate functions in serial.
    #       This will happen on Android and Windows systems
    import MaxwellNRPy_py_dir.MaxwellVacuum_C_kernels_codegen as MCk
    import grid as gri
    for param in paramslist:
        gri.glb_gridfcs_list = []
        MCk.MaxwellVacuum_C_kernels_codegen_onepart(params=param)
    NRPyEnvVars = [] # Reset NRPyEnvVars in case multiprocessing wrote to it and failed.
# # Steps 3.d.ii-iv, alternate: As fallback, evaluate functions in serial.
# #       This will happen on Android and Windows systems

print("(BENCH) Finished C kernel codegen for MaxwellVacuum in "+str(time.time()-start)+" seconds.\n")


# In[8]:


get_ipython().run_line_magic('tb', '')


# <a id='cclfiles'></a>
# 
# # Step 3: ETK `ccl` file generation \[Back to [top](#toc)\]
# $$\label{cclfiles}$$
# 
# The Einstein Toolkit (ETK) ccl files contain runtime parameters (`param.ccl`), registered gridfunctions (`interface.ccl`), and function scheduling (`schedule.ccl`). As parameters and gridfunctions are registered with NRPy+ when the C-code kernels are generated, and this generation occurs on separate processes in parallel, we store the entire NRPy+ environment for *each* process. This results in a tremendous amount of duplication, which is sorted out next. Once all duplicated environment variables (e.g., registered gridfunctions) are removed, we replace the current NRPy+ environment with the new one, by setting `gri.glb_gridfcs_list[],par.glb_params_list[],par.glb_Cparams_list[]`.

# In[9]:


# Store all NRPy+ environment variables to file so NRPy+ environment from within this subprocess can be easily restored
import pickle                   # Standard Python module for converting arbitrary data structures to a uniform format.
import grid as gri              # NRPy+: Functions having to do with numerical grids
import NRPy_param_funcs as par  # NRPy+: Parameter interface

if len(NRPyEnvVars) > 0:
    # https://www.pythonforthelab.com/blog/storing-binary-data-and-serializing/
    grfcs_list = []
    param_list = []
    Cparm_list = []

    for WhichParamSet in NRPyEnvVars[0]:
        # gridfunctions
        i=0
        # print("Length of WhichParamSet:",str(len(WhichParamSet)))
        num_elements = pickle.loads(WhichParamSet[i]); i+=1
        for lst in range(num_elements):
            grfcs_list.append(gri.glb_gridfc(gftype=pickle.loads(WhichParamSet[i+0]),
                                             name  =pickle.loads(WhichParamSet[i+1]),
                                             rank  =pickle.loads(WhichParamSet[i+2]),
                                             DIM   =pickle.loads(WhichParamSet[i+3]))) ; i+=4
        # parameters
        num_elements = pickle.loads(WhichParamSet[i]); i+=1
        for lst in range(num_elements):
            param_list.append(par.glb_param(type      =pickle.loads(WhichParamSet[i+0]),
                                            module    =pickle.loads(WhichParamSet[i+1]),
                                            parname   =pickle.loads(WhichParamSet[i+2]),
                                            defaultval=pickle.loads(WhichParamSet[i+3]))) ; i+=4
        # Cparameters
        num_elements = pickle.loads(WhichParamSet[i]); i+=1
        for lst in range(num_elements):
            Cparm_list.append(par.glb_Cparam(type      =pickle.loads(WhichParamSet[i+0]),
                                             module    =pickle.loads(WhichParamSet[i+1]),
                                             parname   =pickle.loads(WhichParamSet[i+2]),
                                             defaultval=pickle.loads(WhichParamSet[i+3]))) ; i+=4
    grfcs_list_uniq = []
    for gf_ntuple_stored in grfcs_list:
        found_gf = False
        for gf_ntuple_new in grfcs_list_uniq:
            if gf_ntuple_new == gf_ntuple_stored:
                found_gf = True
        if found_gf == False:
            grfcs_list_uniq.append(gf_ntuple_stored)

    param_list_uniq = []
    for pr_ntuple_stored in param_list:
        found_pr = False
        for pr_ntuple_new in param_list_uniq:
            if pr_ntuple_new == pr_ntuple_stored:
                found_pr = True
        if found_pr == False:
            param_list_uniq.append(pr_ntuple_stored)

    # Set glb_paramsvals_list:
    # Step 1: Reset all paramsvals to their defaults
    par.glb_paramsvals_list = []
    for parm in param_list_uniq:
        par.glb_paramsvals_list.append(parm.defaultval)

    Cparm_list_uniq = []
    for Cp_ntuple_stored in Cparm_list:
        found_Cp = False
        for Cp_ntuple_new in Cparm_list_uniq:
            if Cp_ntuple_new == Cp_ntuple_stored:
                found_Cp = True
        if found_Cp == False:
            Cparm_list_uniq.append(Cp_ntuple_stored)

    gri.glb_gridfcs_list = []
    par.glb_params_list  = []
    par.glb_Cparams_list = []

    gri.glb_gridfcs_list = grfcs_list_uniq
    par.glb_params_list  = param_list_uniq
    par.glb_Cparams_list = Cparm_list_uniq


# <a id='paramccl'></a>
# 
# ## Step 3.a: `param.ccl`: specify free parameters within `MaxwellVacuum` \[Back to [top](#toc)\]
# $$\label{paramccl}$$
# 
# All parameters necessary to evolve the right-hand side (RHS) expressions of Maxwell's equations are registered within NRPy+; we use this information to automatically generate `param.ccl`. NRPy+ also specifies default values for each parameter. 
# 
# More information on `param.ccl` syntax can be found in the [official Einstein Toolkit documentation](http://einsteintoolkit.org/usersguide/UsersGuide.html#x1-184000D2.3).

# In[10]:


get_ipython().run_cell_magic('writefile', '$MaxwellVacuumdir/MaxwellVacuum_ETK_ccl_files_codegen.py', '\n# Step 1: Import needed core NRPy+ modules\nimport NRPy_param_funcs as par # NRPy+: Parameter interface\nimport grid as gri             # NRPy+: Functions having to do with numerical grids\nimport os, sys                 # Standard Python modules for multiplatform OS-level functions\n\ndef keep_param__return_type(paramtuple):\n    keep_param = True # We\'ll not set some parameters in param.ccl;\n                      #   e.g., those that should be #define\'d like M_PI.\n    typestring = ""\n    # Separate thorns within the ETK take care of grid/coordinate parameters;\n    #   thus we ignore NRPy+ grid/coordinate parameters:\n    if paramtuple.module == "grid" or paramtuple.module == "reference_metric":\n        keep_param = False\n\n    partype = paramtuple.type\n    if partype == "bool":\n        typestring += "BOOLEAN "\n    elif partype == "REAL":\n        if paramtuple.defaultval != 1e300: # 1e300 is a magic value indicating that the C parameter should be mutable\n            typestring += "CCTK_REAL "\n        else:\n            keep_param = False\n    elif partype == "int":\n        typestring += "CCTK_INT "\n    elif partype == "#define":\n        keep_param = False\n    elif partype == "char":\n        # FIXME: char/string parameter types should in principle be supported\n        print("Error: parameter "+paramtuple.module+"::"+paramtuple.parname+\n              " has unsupported type: \\""+ paramtuple.type + "\\"")\n        sys.exit(1)\n    else:\n        print("Error: parameter "+paramtuple.module+"::"+paramtuple.parname+\n              " has unsupported type: \\""+ paramtuple.type + "\\"")\n        sys.exit(1)\n    return keep_param, typestring\n\ndef output_param_ccl(ThornName="MaxwellVacuum"):\n    with open(os.path.join(ThornName,"param.ccl"), "w") as file:\n        file.write("""\n# This param.ccl file was automatically generated by NRPy+.\n#   You are advised against modifying it directly; instead\n#   modify the Python code that generates it.\n\nshares: MethodOfLines\n\n#EXTENDS CCTK_KEYWORD evolution_method "evolution_method"\n#{\n#  "MaxwellVacuum" :: ""\n#}\n\nrestricted:\n\nCCTK_INT FD_order "Finite-differencing order"\n{\\n""")\n        FDorders = []\n        for _root, _dirs, files in os.walk(os.path.join(ThornName,"src")): # _root,_dirs unused.\n            for Cfilename in files:\n                if (".h" in Cfilename) and ("RHSs" in Cfilename) and ("intrinsics" not in Cfilename):\n                    array = Cfilename.replace(".","_").split("_")\n                    FDorders.append(int(array[-2]))\n        FDorders.sort()\n        for order in FDorders:\n            file.write(" "+str(order)+":"+str(order)+"   :: \\"finite-differencing order = "+str(order)+"\\"\\n")\n        FDorder_default = 4\n        if FDorder_default not in FDorders:\n            print("WARNING: 4th-order FD kernel was not output!?! Changing default FD order to "+str(FDorders[0]))\n            FDorder_default = FDorders[0]\n        file.write("} "+str(FDorder_default)+ "\\n\\n") # choose 4th order by default, consistent with ML_Maxwell\n        paramccl_str = ""\n        for i in range(len(par.glb_Cparams_list)):\n            # keep_param is a boolean indicating whether we should accept or reject\n            #    the parameter. singleparstring will contain the string indicating\n            #    the variable type.\n            keep_param, singleparstring = keep_param__return_type(par.glb_Cparams_list[i])\n\n            if keep_param:\n                parname = par.glb_Cparams_list[i].parname\n                partype = par.glb_Cparams_list[i].type\n                singleparstring += parname + " \\""+ parname +" (see NRPy+ for parameter definition)\\"\\n"\n                singleparstring += "{\\n"\n                if partype != "bool":\n                    singleparstring += " *:* :: \\"All values accepted. NRPy+ does not restrict the allowed ranges of parameters yet.\\"\\n"\n                singleparstring += "} "+str(par.glb_Cparams_list[i].defaultval)+"\\n\\n"\n\n                paramccl_str += singleparstring\n        file.write(paramccl_str)\n')


# <a id='interfaceccl'></a>
# 
# ## Step 3.b: `interface.ccl`: define needed gridfunctions; provide keywords denoting what this thorn provides and what it should inherit from other thorns \[Back to [top](#toc)\]
# $$\label{interfaceccl}$$
# 
# `interface.ccl` declares all gridfunctions and determines how `MaxwellVacuum` interacts with other Einstein Toolkit thorns.
# 
# The [official Einstein Toolkit (Cactus) documentation](http://einsteintoolkit.org/usersguide/UsersGuide.html#x1-179000D2.2) defines what must/should be included in an `interface.ccl` file.

# In[11]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_ETK_ccl_files_codegen.py', '\n# First construct lists of the basic gridfunctions used in NRPy+.\n#    Each type will be its own group in MaxwellVacuum.\nevol_gfs_list    = []\nauxevol_gfs_list = []\naux_gfs_list     = []\nfor i in range(len(gri.glb_gridfcs_list)):\n    if gri.glb_gridfcs_list[i].gftype == "EVOL":\n        evol_gfs_list.append(   gri.glb_gridfcs_list[i].name+"GF")\n\n    if gri.glb_gridfcs_list[i].gftype == "AUX":\n        aux_gfs_list.append(    gri.glb_gridfcs_list[i].name+"GF")\n\n    if gri.glb_gridfcs_list[i].gftype == "AUXEVOL":\n        auxevol_gfs_list.append(gri.glb_gridfcs_list[i].name+"GF")\n\n# NRPy+\'s finite-difference code generator assumes gridfunctions\n#    are alphabetized; not sorting may result in unnecessary\n#    cache misses.\nevol_gfs_list.sort()\naux_gfs_list.sort()\nauxevol_gfs_list.sort()\nrhs_list = []\nfor gf in evol_gfs_list:\n    rhs_list.append(gf.replace("GF","")+"_rhsGF")\n\ndef output_interface_ccl(ThornName="MaxwellVacuum"):\n    outstr = """\n# This interface.ccl file was automatically generated by NRPy+.\n#   You are advised against modifying it directly; instead\n#   modify the Python code that generates it.\n\n# With "implements", we give our thorn its unique name.\nimplements: MaxwellVacuum\n\n# By "inheriting" other thorns, we tell the Toolkit that we\n#   will rely on variables/function that exist within those\n#   functions.\ninherits: Boundary grid MethodofLines\\n"""\n\n    outstr += """\n# Needed functions and #include\'s:\nUSES INCLUDE: Symmetry.h\nUSES INCLUDE: Boundary.h\n\n# Needed Method of Lines function\nCCTK_INT FUNCTION MoLRegisterEvolvedGroup(CCTK_INT IN EvolvedIndex, \\\n                                          CCTK_INT IN RHSIndex)\nREQUIRES FUNCTION MoLRegisterEvolvedGroup\n\n# Needed Boundary Conditions function\nCCTK_INT FUNCTION GetBoundarySpecification(CCTK_INT IN size, CCTK_INT OUT ARRAY nboundaryzones, CCTK_INT OUT ARRAY is_internal, CCTK_INT OUT ARRAY is_staggered, CCTK_INT OUT ARRAY shiftout)\nUSES FUNCTION GetBoundarySpecification\n\nCCTK_INT FUNCTION SymmetryTableHandleForGrid(CCTK_POINTER_TO_CONST IN cctkGH)\nUSES FUNCTION SymmetryTableHandleForGrid\n\nCCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, CCTK_STRING IN var_name, CCTK_STRING IN bc_name)\nUSES FUNCTION Boundary_SelectVarForBC\n\n# Needed for EinsteinEvolve/NewRad outer boundary condition driver:\nCCTK_INT FUNCTION                         \\\\\n    NewRad_Apply                          \\\\\n        (CCTK_POINTER_TO_CONST IN cctkGH, \\\\\n         CCTK_REAL ARRAY IN var,          \\\\\n         CCTK_REAL ARRAY INOUT rhs,       \\\\\n         CCTK_REAL IN var0,               \\\\\n         CCTK_REAL IN v0,                 \\\\\n         CCTK_INT IN radpower)\nREQUIRES FUNCTION NewRad_Apply\n\n# Tell the Toolkit that we want all gridfunctions\n#    to be visible to other thorns by using\n#    the keyword "public". Note that declaring these\n#    gridfunctions *does not* allocate memory for them;\n#    that is done by the schedule.ccl file.\n\npublic:\n"""\n\n    # Next we declare gridfunctions based on their corresponding gridfunction groups as registered within NRPy+\n\n    def output_list_of_gfs(gfs_list,description="User did not provide description"):\n        gfsstr = "    "\n        for i in range(len(gfs_list)):\n            gfsstr += gfs_list[i]\n            if i != len(gfs_list)-1:\n                gfsstr += "," # This is a comma-separated list of gridfunctions\n            else:\n                gfsstr += "\\n} \\""+description+"\\"\\n\\n"\n        return gfsstr\n    # First EVOL type:\n    outstr += "CCTK_REAL evol_variables type = GF Timelevels=3\\n{\\n"\n    outstr += output_list_of_gfs(evol_gfs_list,"Maxwell evolved gridfunctions")\n    # Second EVOL right-hand-sides\n    outstr += "CCTK_REAL evol_variables_rhs type = GF Timelevels=1 TAGS=\\\'InterpNumTimelevels=1 prolongation=\\"none\\"\\\'\\n{\\n"\n    outstr += output_list_of_gfs(rhs_list,"right-hand-side storage for Maxwell evolved gridfunctions")\n    # Then AUX type:\n    outstr += "CCTK_REAL aux_variables type = GF Timelevels=3\\n{\\n"\n    outstr += output_list_of_gfs(aux_gfs_list,"Auxiliary gridfunctions for Maxwell diagnostics")\n    # Finally, AUXEVOL type:\n#     outstr += "CCTK_REAL auxevol_variables type = GF Timelevels=1 TAGS=\\\'InterpNumTimelevels=1 prolongation=\\"none\\"\\\'\\n{\\n"\n#     outstr += output_list_of_gfs(auxevol_gfs_list,"Auxiliary gridfunctions needed for evaluating the Maxwell RHSs")\n\n    with open(os.path.join(ThornName,"interface.ccl"), "w") as file:\n        file.write(outstr.replace("MaxwellVacuum",ThornName))\n')


# <a id='scheduleccl'></a>
# 
# ## Step 3.c: `schedule.ccl`: schedule all functions used within `MaxwellVacuum`, specify data dependencies within said functions, and allocate memory for gridfunctions \[Back to [top](#toc)\]
# $$\label{scheduleccl}$$
# 
# Official documentation on constructing ETK `schedule.ccl` files is found [here](http://einsteintoolkit.org/usersguide/UsersGuide.html#x1-187000D2.4).

# In[12]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_ETK_ccl_files_codegen.py', '\ndef output_schedule_ccl(ThornName="MaxwellVacuum"):\n    outstr = """\n# This schedule.ccl file was automatically generated by NRPy+.\n#   You are advised against modifying it directly; instead\n#   modify the Python code that generates it.\n\n# Next allocate storage for all 3 gridfunction groups used in MaxwellVacuum\nSTORAGE: evol_variables[3]     # Evolution variables\nSTORAGE: evol_variables_rhs[1] # Variables storing right-hand-sides\nSTORAGE: aux_variables[3]      # Diagnostics variables\n\n# The following scheduler is based on Lean/LeanMaxwellMoL/schedule.ccl\n\nschedule MaxwellVacuum_Banner at STARTUP\n{\n  LANG: C\n  OPTIONS: meta\n} "Output ASCII art banner"\n\nschedule MaxwellVacuum_Symmetry_registration at BASEGRID\n{\n  LANG: C\n  OPTIONS: Global\n} "Register symmetries, the CartGrid3D way."\n\nschedule MaxwellVacuum_zero_rhss at BASEGRID after MaxwellVacuum_Symmetry_registration\n{\n  LANG: C\n} "Idea from Lean: set all rhs functions to zero to prevent spurious nans"\n\n# MoL: registration\n\nschedule MaxwellVacuum_MoL_registration in MoL_Register\n{\n  LANG: C\n  OPTIONS: META\n} "Register variables for MoL"\n\n# MoL: compute RHSs, etc\n\nschedule MaxwellVacuum_RHSs in MoL_CalcRHS as MaxwellVacuum_RHS\n{\n  LANG: C\n} "MoL: Evaluate Maxwell RHSs"\n\nschedule MaxwellVacuum_NewRad in MoL_CalcRHS after MaxwellVacuum_RHS\n{\n  LANG: C\n} "NewRad boundary conditions, scheduled right after RHS eval."\n\nschedule MaxwellVacuum_BoundaryConditions_evolved_gfs in MoL_PostStep\n{\n  LANG: C\n  OPTIONS: LEVEL\n  SYNC: evol_variables\n} "Apply boundary conditions and perform AMR+interprocessor synchronization"\n\nschedule GROUP ApplyBCs as MaxwellVacuum_ApplyBCs in MoL_PostStep after MaxwellVacuum_BoundaryConditions_evolved_gfs\n{\n} "Group for applying boundary conditions"\n\n# Compute divergence and Gamma constraints\n\nschedule MaxwellVacuum_constraints in MoL_PseudoEvolution\n{\n  LANG: C\n  OPTIONS: Local\n} "Compute Maxwell (divergence and Gamma) constraints"\n"""\n    with open(os.path.join(ThornName,"schedule.ccl"), "w") as file:\n        file.write(outstr.replace("MaxwellVacuum",ThornName))\n')


# In[13]:


import MaxwellNRPy_py_dir.MaxwellVacuum_ETK_ccl_files_codegen as cclgen

ThornName="MaxwellVacuum"
cclgen.output_param_ccl(ThornName)
cclgen.output_interface_ccl(ThornName)
cclgen.output_schedule_ccl(ThornName)


# <a id='cdrivers'></a>
# 
# # Step 4: C driver functions for ETK registration & NRPy+ generated kernels \[Back to [top](#toc)\]
# $$\label{cdrivers}$$
# 
# Now that we have constructed the basic C code kernels and the needed Einstein Toolkit `ccl` files, we next write the driver functions for registering `MaxwellVacuum` within the Toolkit and the C code kernels. Each of these driver functions will be called directly from the thorn's [`schedule.ccl`](#scheduleccl) in the ETK.
# 
# <a id='etkfunctions'></a>
# ## Step 4.a: Needed ETK functions: Banner, Symmetry registration, Parameter sanity check,  Method of Lines (`MoL`) registration, Boundary condition \[Back to [top](#toc)\]
# $$\label{etkfunctions}$$
# 
# ### To-do: Parameter sanity check function. E.g., error should be thrown if `cctk_nghostzones[]` is set too small for the chosen finite-differencing order within NRPy+.

# In[14]:


get_ipython().run_cell_magic('writefile', '$MaxwellVacuumdir/MaxwellVacuum_C_drivers_codegen.py', '\n# Step 1: Import needed core NRPy+ and Python modules\nfrom outputC import lhrh         # NRPy+: Core C code output module\nimport NRPy_param_funcs as par   # NRPy+: Parameter interface\nimport finite_difference as fin  # NRPy+: Finite difference C code generation module\nimport grid as gri               # NRPy+: Functions having to do with numerical grids\nimport loop as lp                # NRPy+: Generate C code loops\nimport indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\nimport os, sys                   # Standard Python modules for multiplatform OS-level functions\n# We need the function keep_param__return_type() from this module:\nimport MaxwellNRPy_py_dir.MaxwellVacuum_ETK_ccl_files_codegen as ccl\n\nmake_code_defn_list = []\ndef append_to_make_code_defn_list(filename):\n    if filename not in make_code_defn_list:\n        make_code_defn_list.append(filename)\n    return filename\n')


# In[15]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_drivers_codegen.py', '\ndef driver_C_codes(Csrcdict, ThornName,\n                   rhs_list,evol_gfs_list,aux_gfs_list,auxevol_gfs_list):\n    # First the ETK banner code, proudly showing the NRPy+ banner\n    import NRPy_logo as logo\n    outstr = """\n#include <stdio.h>\n\nvoid MaxwellVacuum_Banner()\n{\n    """\n    logostr = logo.print_logo(print_to_stdout=False)\n    outstr += "printf(\\"MaxwellVacuum: another Einstein Toolkit thorn generated by\\\\n\\");\\n"\n    for line in logostr.splitlines():\n        outstr += "    printf(\\""+line+"\\\\n\\");\\n"\n    outstr += "}\\n"\n\n    # Finally add C code string to dictionaries (Python dictionaries are immutable)\n    # Add C code string to dictionary (Python dictionaries are immutable)\n    Csrcdict[append_to_make_code_defn_list("Banner.c")] = outstr.replace("MaxwellVacuum",ThornName)\n')


# In[16]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_drivers_codegen.py', '\n    # Next MaxwellVacuum_Symmetry_registration(): Register symmetries\n\n    full_gfs_list = []\n    full_gfs_list.extend(evol_gfs_list)\n    full_gfs_list.extend(auxevol_gfs_list)\n    full_gfs_list.extend(aux_gfs_list)\n\n    outstr = """\n#include "cctk.h"\n#include "cctk_Arguments.h"\n#include "cctk_Parameters.h"\n#include "Symmetry.h"\n\nvoid MaxwellVacuum_Symmetry_registration(CCTK_ARGUMENTS)\n{\n  DECLARE_CCTK_ARGUMENTS;\n  DECLARE_CCTK_PARAMETERS;\n\n  // Stores gridfunction parity across x=0, y=0, and z=0 planes, respectively\n  int sym[3];\n\n  // Next register parities for each gridfunction based on its name\n  //    (to ensure this algorithm is robust, gridfunctions with integers\n  //     in their base names are forbidden in NRPy+).\n"""\n    outstr += ""\n    for gfname in full_gfs_list:\n        gfname_without_GFsuffix = gfname[:-2]\n        outstr += """\n      // Default to scalar symmetry:\n      sym[0] = 1; sym[1] = 1; sym[2] = 1;\n      // Now modify sym[0], sym[1], and/or sym[2] as needed\n      //    to account for gridfunction parity across\n      //    x=0, y=0, and/or z=0 planes, respectively\n"""\n        # If gridfunction name does not end in a digit, by NRPy+ syntax, it must be a scalar\n        if gfname_without_GFsuffix[len(gfname_without_GFsuffix) - 1].isdigit() == False:\n            outstr += "      // (this gridfunction is a scalar -- no need to change default sym[]\'s!)\\n"\n        elif len(gfname_without_GFsuffix) > 2:\n            # Rank-1 indexed expression (e.g., vector)\n            if gfname_without_GFsuffix[len(gfname_without_GFsuffix) - 2].isdigit() == False:\n                if int(gfname_without_GFsuffix[-1]) > 2:\n                    print("Error: Found invalid gridfunction name: "+gfname)\n                    sys.exit(1)\n                symidx = gfname_without_GFsuffix[-1]\n                if int(symidx) < 3:   outstr += "      sym[" + symidx + "] = -1;\\n"\n            # Rank-2 indexed expression\n            elif gfname_without_GFsuffix[len(gfname_without_GFsuffix) - 2].isdigit() == True:\n                if len(gfname_without_GFsuffix) > 3 and gfname_without_GFsuffix[len(gfname_without_GFsuffix) - 3].isdigit() == True:\n                    print("Error: Found a Rank-3 or above gridfunction: "+gfname+", which is at the moment unsupported.")\n                    print("It should be easy to support this if desired.")\n                    sys.exit(1)\n                symidx0 = gfname_without_GFsuffix[-2]\n                if int(symidx0) >= 0: outstr += "      sym[" + symidx0 + "] *= -1;\\n"\n                symidx1 = gfname_without_GFsuffix[-1]\n                if int(symidx1) >= 0: outstr += "      sym[" + symidx1 + "] *= -1;\\n"\n        else:\n            print("Don\'t know how you got this far with a gridfunction named "+gfname+", but I\'ll take no more of this nonsense.")\n            print("   Please follow best-practices and rename your gridfunction to be more descriptive")\n            sys.exit(1)\n        outstr += "      SetCartSymVN(cctkGH, sym, \\"MaxwellVacuum::" + gfname + "\\");\\n"\n    outstr += "}\\n"\n\n    # Add C code string to dictionary (Python dictionaries are immutable)\n    Csrcdict[append_to_make_code_defn_list("Symmetry_registration_oldCartGrid3D.c")] = \\\n        outstr.replace("MaxwellVacuum",ThornName)\n')


# In[17]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_drivers_codegen.py', '\n    # Next set RHSs to zero\n    outstr = """\n#include "cctk.h"\n#include "cctk_Arguments.h"\n#include "cctk_Parameters.h"\n#include "Symmetry.h"\n\nvoid MaxwellVacuum_zero_rhss(CCTK_ARGUMENTS)\n{\n  DECLARE_CCTK_ARGUMENTS;\n  DECLARE_CCTK_PARAMETERS;\n"""\n    set_rhss_to_zero = ""\n    for gf in rhs_list:\n        set_rhss_to_zero += gf+"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)] = 0.0;\\n"\n\n    outstr += lp.loop(["i2","i1","i0"],["0", "0", "0"],\n                      ["cctk_lsh[2]","cctk_lsh[1]","cctk_lsh[0]"],\n                      ["1","1","1"],\n                      ["#pragma omp parallel for","","",],"",set_rhss_to_zero)\n    outstr += "}\\n"\n    # Add C code string to dictionary (Python dictionaries are immutable)\n    Csrcdict[append_to_make_code_defn_list("zero_rhss.c")] = outstr.replace("MaxwellVacuum",ThornName)\n')


# In[18]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_drivers_codegen.py', '\n    # Next registration with the Method of Lines thorn\n    outstr = """\n//--------------------------------------------------------------------------\n// Register with the Method of Lines time stepper\n// (MoL thorn, found in arrangements/CactusBase/MoL)\n// MoL documentation located in arrangements/CactusBase/MoL/doc\n//--------------------------------------------------------------------------\n#include <stdio.h>\n\n#include "cctk.h"\n#include "cctk_Parameters.h"\n#include "cctk_Arguments.h"\n\n#include "Symmetry.h"\n\nvoid MaxwellVacuum_MoL_registration(CCTK_ARGUMENTS)\n{\n  DECLARE_CCTK_ARGUMENTS;\n  DECLARE_CCTK_PARAMETERS;\n\n  CCTK_INT ierr = 0, group, rhs;\n\n  // Register evolution & RHS gridfunction groups with MoL, so it knows\n\n  group = CCTK_GroupIndex("MaxwellVacuum::evol_variables");\n  rhs = CCTK_GroupIndex("MaxwellVacuum::evol_variables_rhs");\n  ierr += MoLRegisterEvolvedGroup(group, rhs);\n\n  if (ierr) CCTK_ERROR("Problems registering with MoL");\n}\n"""\n    # Add C code string to dictionary (Python dictionaries are immutable)\n    Csrcdict[append_to_make_code_defn_list("MoL_registration.c")] = outstr.replace("MaxwellVacuum",ThornName)\n')


# In[19]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_drivers_codegen.py', '\n    # Next register with the boundary conditions thorns.\n    # PART 1: Set BC type to "none" for all variables\n    # Since we choose NewRad boundary conditions, we must register all\n    #   gridfunctions to have boundary type "none". This is because\n    #   NewRad is seen by the rest of the Toolkit as a modification to the\n    #   RHSs.\n\n    # This code is based on Kranc\'s McLachlan/ML_Maxwell/src/Boundaries.cc code.\n    outstr = """\n#include "cctk.h"\n#include "cctk_Arguments.h"\n#include "cctk_Parameters.h"\n#include "cctk_Faces.h"\n#include "util_Table.h"\n#include "Symmetry.h"\n\n// Set `none` boundary conditions on Maxwell RHSs, as these are set via NewRad.\nvoid MaxwellVacuum_BoundaryConditions_evolved_gfs(CCTK_ARGUMENTS)\n{\n  DECLARE_CCTK_ARGUMENTS;\n  DECLARE_CCTK_PARAMETERS;\n\n  CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0;\n"""\n    for gf in evol_gfs_list:\n        outstr += """\n  ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, "MaxwellVacuum::"""+gf+"""", "none");\n  if (ierr < 0) CCTK_ERROR("Failed to register BC for MaxwellVacuum::"""+gf+"""!");\n"""\n    outstr += """\n}\n\n// Set `none` boundary conditions on Maxwell constraints\nvoid MaxwellVacuum_BoundaryConditions_aux_gfs(CCTK_ARGUMENTS) {\n  DECLARE_CCTK_ARGUMENTS;\n  DECLARE_CCTK_PARAMETERS;\n\n  CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0;\n\n"""\n    for gf in aux_gfs_list:\n        outstr += """\n  ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, cctk_nghostzones[0], -1, "MaxwellVacuum::"""+gf+"""", "none");\n  if (ierr < 0) CCTK_ERROR("Failed to register BC for MaxwellVacuum::"""+gf+"""!");\n"""\n    outstr += "}\\n"\n\n    # Add C code string to dictionary (Python dictionaries are immutable)\n    Csrcdict[append_to_make_code_defn_list("BoundaryConditions.c")] = outstr.replace("MaxwellVacuum",ThornName)\n\n    # PART 2: Set C code for calling NewRad BCs\n    #   As explained in lean_public/LeanMaxwellMoL/src/calc_mwev_rhs.F90,\n    #   the function NewRad_Apply takes the following arguments:\n    #   NewRad_Apply(cctkGH, var, rhs, var0, v0, radpower),\n    #     which implement the boundary condition:\n    #       var  =  var_at_infinite_r + u(r-var_char_speed*t)/r^var_radpower\n    #  Obviously for var_radpower>0, var_at_infinite_r is the value of\n    #    the variable at r->infinity. var_char_speed is the propagation\n    #    speed at the outer boundary, and var_radpower is the radial\n    #    falloff rate.\n\n    outstr = """\n#include <math.h>\n\n#include "cctk.h"\n#include "cctk_Arguments.h"\n#include "cctk_Parameters.h"\n\nvoid MaxwellVacuum_NewRad(CCTK_ARGUMENTS) {\n  DECLARE_CCTK_ARGUMENTS;\n  DECLARE_CCTK_PARAMETERS;\n\n"""\n    for gf in evol_gfs_list:\n        var_at_infinite_r = "0.0"\n        var_char_speed    = "1.0"\n        var_radpower      = "3.0"\n\n        outstr += "  NewRad_Apply(cctkGH, "+gf+", "+gf.replace("GF","")+"_rhsGF, "+var_at_infinite_r+", "+var_char_speed+", "+var_radpower+");\\n"\n    outstr += "}\\n"\n\n    # Add C code string to dictionary (Python dictionaries are immutable)\n    Csrcdict[append_to_make_code_defn_list("BoundaryCondition_NewRad.c")] = outstr.replace("MaxwellVacuum",ThornName)\n')


# <a id='mwevrhss'></a>
# 
# ## Step 4.b: Evaluate Maxwell right-hand-sides (RHSs) \[Back to [top](#toc)\]
# $$\label{mwevrhss}$$

# In[20]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_drivers_codegen.py', '\n    ###########################\n    ###########################\n    # Maxwell_RHSs\n    ###########################\n    common_includes = """\n#include <math.h>\n#include "cctk.h"\n#include "cctk_Arguments.h"\n#include "cctk_Parameters.h"\n#include "SIMD/SIMD_intrinsics.h"\n"""\n    common_preloop = """\n    DECLARE_CCTK_ARGUMENTS;\n    const CCTK_REAL NOSIMDinvdx0 = 1.0/CCTK_DELTA_SPACE(0);\n    const REAL_SIMD_ARRAY invdx0 = ConstSIMD(NOSIMDinvdx0);\n    const CCTK_REAL NOSIMDinvdx1 = 1.0/CCTK_DELTA_SPACE(1);\n    const REAL_SIMD_ARRAY invdx1 = ConstSIMD(NOSIMDinvdx1);\n    const CCTK_REAL NOSIMDinvdx2 = 1.0/CCTK_DELTA_SPACE(2);\n    const REAL_SIMD_ARRAY invdx2 = ConstSIMD(NOSIMDinvdx2);\n"""\n')


# In[21]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_drivers_codegen.py', '\n    path = os.path.join(ThornName,"src")\n    MaxwellVacuum_src_filelist = []\n    for _root, _dirs, files in os.walk(path):  # _root, _dirs unused.\n        for filename in files:\n            MaxwellVacuum_src_filelist.append(filename)\n    MaxwellVacuum_src_filelist.sort() # Sort the list in place.\n\n    Maxwell_FD_orders_output = []\n    for filename in MaxwellVacuum_src_filelist:\n        if "Maxwell_RHSs_" in filename:\n            array = filename.replace(".","_").split("_")\n            FDorder =  int(array[-2])\n            if FDorder not in Maxwell_FD_orders_output:\n                Maxwell_FD_orders_output.append(FDorder)\n    Maxwell_FD_orders_output.sort()\n\n    ###########################\n    # Output Maxwell RHSs driver function\n    outstr = common_includes\n    for filename in MaxwellVacuum_src_filelist:\n        if ("Maxwell_RHSs_" in filename) and (".h" in filename):\n            outstr += """extern void """ + ThornName+"_"+filename.replace(".h", "(CCTK_ARGUMENTS);") + "\\n"\n\n    outstr += """\nvoid MaxwellVacuum_RHSs(CCTK_ARGUMENTS) {\n    DECLARE_CCTK_ARGUMENTS;\n\n    const CCTK_INT *FD_order = CCTK_ParameterGet("FD_order","MaxwellVacuum",NULL);\n\n"""\n    for filename in MaxwellVacuum_src_filelist:\n        if ("Maxwell_RHSs_" in filename) and (".h" in filename):\n            array = filename.replace(".", "_").split("_")\n            outstr += "    if(*FD_order == " + str(array[-2]) + ") {\\n"\n            outstr += "        " + ThornName+"_"+filename.replace(".h", "(CCTK_PASS_CTOC);") + "\\n"\n            outstr += "    }\\n"\n    outstr += "} // END FUNCTION\\n"\n    # Add C code string to dictionary (Python dictionaries are immutable)\n    Csrcdict[append_to_make_code_defn_list("Maxwell_RHSs.c")] = outstr.replace("MaxwellVacuum", ThornName)\n\n    def SIMD_declare_C_params():\n        SIMD_declare_C_params_str = ""\n        for i in range(len(par.glb_Cparams_list)):\n            # keep_param is a boolean indicating whether we should accept or reject\n            #    the parameter. singleparstring will contain the string indicating\n            #    the variable type.\n            keep_param, singleparstring = ccl.keep_param__return_type(par.glb_Cparams_list[i])\n\n            if (keep_param) and ("CCTK_REAL" in singleparstring):\n                parname = par.glb_Cparams_list[i].parname\n                SIMD_declare_C_params_str += "    const "+singleparstring + "*NOSIMD"+parname+\\\n                                             " = CCTK_ParameterGet(\\""+parname+"\\",\\"MaxwellVacuum\\",NULL);\\n"\n                SIMD_declare_C_params_str += "    const REAL_SIMD_ARRAY "+parname+" = ConstSIMD(*NOSIMD"+parname+");\\n"\n        return SIMD_declare_C_params_str\n\n    # Create functions for the largest C kernels (Maxwell RHSs and Ricci) and output\n    #    the .h files to .c files with function wrappers; delete original .h files\n    path = os.path.join(ThornName, "src")\n    for filename in MaxwellVacuum_src_filelist:\n        if ("Maxwell_RHSs_" in filename) and (".h" in filename):\n            outstr = common_includes + "void MaxwellVacuum_"+filename.replace(".h","")+"(CCTK_ARGUMENTS) {\\n"\n            outstr += common_preloop+SIMD_declare_C_params()\n            with open(os.path.join(path,filename), "r") as currfile:\n                outstr += currfile.read()\n            # Now that we\'ve inserted the contents of the kernel into this file,\n            #     we delete the file containing the kernel\n            os.remove(os.path.join(path,filename))\n            outstr += "} // END FUNCTION\\n"\n            # Add C code string to dictionary (Python dictionaries are immutable)\n            Csrcdict[append_to_make_code_defn_list(filename.replace(".h",".c"))] = outstr.replace("MaxwellVacuum",ThornName)\n')


# <a id='diagnostics'></a>
# 
# ## Step 4.c: Diagnostics: Computing the divergence constraint \[Back to [top](#toc)\]
# $$\label{diagnostics}$$
# 
# The divergence constraint is a useful diagnostics of a calculation's health. Here we construct the driver function.

# In[22]:


get_ipython().run_cell_magic('writefile', '-a $MaxwellVacuumdir/MaxwellVacuum_C_drivers_codegen.py', '\n    # Next, the driver for computing the Maxwell Hamiltonian & momentum constraints\n    outstr = """\n#include <math.h>\n\n#include "cctk.h"\n#include "cctk_Arguments.h"\n#include "cctk_Parameters.h"\n\nvoid MaxwellVacuum_constraints(CCTK_ARGUMENTS) {\n    DECLARE_CCTK_ARGUMENTS;\n    DECLARE_CCTK_PARAMETERS;\n\n    const CCTK_REAL invdx0 = 1.0/CCTK_DELTA_SPACE(0);\n    const CCTK_REAL invdx1 = 1.0/CCTK_DELTA_SPACE(1);\n    const CCTK_REAL invdx2 = 1.0/CCTK_DELTA_SPACE(2);\n"""\n    for filename in MaxwellVacuum_src_filelist:\n        if "Maxwell_constraints_" in filename:\n            array = filename.replace(".","_").split("_")\n            outstr += "    if(FD_order == "+str(array[-2])+") {\\n"\n            outstr += "        #include \\""+filename+"\\"\\n"\n            outstr += "    }\\n"\n    outstr += "}\\n"\n\n    # Add C code string to dictionary (Python dictionaries are immutable)\n    Csrcdict[append_to_make_code_defn_list("driver_constraints.c")] = outstr.replace("MaxwellVacuum",ThornName)\n')


# <a id='outcdrivers'></a>
# 
# ## Step 4.d: Output all C driver functions needed for `MaxwellVacuum` \[Back to [top](#toc)\]
# $$\label{outcdrivers}$$
# 
# First we call the above functions (output above to the `MaxwellVacuum_validate.MaxwellVacuum_C_drivers_codegen` Python module) to store all needed driver C files to a Python dictionary, then we simply outputs the dictionary to the appropriate files.

# In[23]:


import MaxwellNRPy_py_dir.MaxwellVacuum_C_drivers_codegen as driver

# The following Python dictionaries consist of a key, which is the filename
#    in the thorn's src/ directory (e.g., "driver_Maxwell_constraints.c"),
#    and a value, which is the corresponding source code, stored as a
#    Python string.
Vac_Csrcdict = {}
Reg_Csrcdict = {}

# We'll need lists of gridfunctions for these driver functions
evol_gfs_list    = cclgen.evol_gfs_list
aux_gfs_list     = cclgen.aux_gfs_list
auxevol_gfs_list = cclgen.auxevol_gfs_list

# Generate driver codes for MaxwellVacuum thorn (i.e., populate the Vac_Csrcdict dictionary)
driver.driver_C_codes(Vac_Csrcdict, "MaxwellVacuum",
                               cclgen.rhs_list,cclgen.evol_gfs_list,cclgen.aux_gfs_list,
                               cclgen.auxevol_gfs_list)

# Next we output the contents of the Reg_Csrcdict and
#   Vac_Csrcdict dictionaries to files in the respective
#   thorns' directories.
for key,val in Vac_Csrcdict.items():
    with open(os.path.join("MaxwellVacuum","src",key),"w") as file:
        file.write(val)


# <a id='makecodedefn'></a>
# 
# ## Step 4.e: `make.code.defn`: List of all C driver functions needed to compile `MaxwellVacuum` \[Back to [top](#toc)\]
# $$\label{makecodedefn}$$
# 
# When constructing each C code driver function above, we called the `append_to_make_code_defn_list()` function, which built a list of each C code driver file. We'll now add each of those files to the `make.code.defn` file, used by the Einstein Toolkit's build system.

# In[24]:


# Finally output the thorns' make.code.defn files, consisting of
#   a list of all C codes in the above dictionaries. This is
#   part of the ETK build system so that these files are output.

def output_make_code_defn(dictionary, ThornName):
    with open(os.path.join(ThornName,"src","make.code.defn"), "w") as file:
        file.write("""
# Main make.code.defn file for thorn """+ThornName+"""

# Source files in this directory
SRCS =""")
        filestring = ""
        list_of_C_driver_files = list(dictionary.keys())
        for i in range(len(list_of_C_driver_files)):
            filestring += "      "+list_of_C_driver_files[i]
            if i != len(list_of_C_driver_files)-1:
                filestring += " \\\n"
            else:
                filestring += "\n"
        file.write(filestring)

output_make_code_defn(Vac_Csrcdict,"MaxwellVacuum")


# <a id='code_validation'></a>
# 
# # Step 5: Code validation \[Back to [top](#toc)\]
# $$\label{code_validation}$$
# 
# Here we will show plots demonstrating good error convergence and proper behavior of error nodes in the systems.

# <a id='convergence'></a>
# 
# ## Step 5.a: Error Convergence \[Back to [top](#toc)\]
# $$\label{convergence}$$
# 
# **Code tests adopting fourth-order finite differencing, coupled to 2nd order Iterative Crank-Nicholson method-of-lines for time integration**
# 
# Inside the directory *`MaxwellVacuum/example_parfiles/`* are the files used for this convergence test:
# **maxwell_toroidaldipole-0.125_OB4.par & maxwell_toroidaldipole-0.0625_OB4.par** : ETK parameter files needed for performing the tests. These parameter files set up a toroidal dipole field propagating in a 3D numerical grid that extends from -4. to +4. along the x-, y-, and z-axes (in units of $c=1$). The parameter files are identical, except the latter has grid resolution that is twice as high as the former (so the errors should drop in the higher resolution case by a factor of $2^2$, since we adopt fourth-order finite differencing coupled to 2nd order Iterative Crank-Nicholson time integration.)
# 
# **Second-order code validation test results:**
# 
# The plot below shows the discrepancy between numerical and exact solutions to x-components of system I $\vec{E}$ and $\vec{A}$ at two different resolutions, at t = 2.0 (to not have errors at the boundary propagate too far inward): dashed is low resolution ($\Delta x_{\rm low}=0.125$) and solid is high resolution ($\Delta x_{\rm high}=0.0625$). Since this test adopts **fourth**-order finite differencing for spatial derivatives and **second**-order Iterative Crank-Nicholson timestepping, we would expect this error to drop by a factor of approximately $(\Delta x_{\rm low}/\Delta x_{\rm high})^2 = (0.125/0.0625)^2 = 2^2=4$ when going from low to high resolution, and after rescaling the error in the low-resolution case, we see that indeed it overlaps the high-resolution result quite nicely, confirming second-order convergence. We note that we also observe convergence for all other evolved variables (in both systems) with a nonzero exact solutions.

# In[25]:


from IPython.display import Image
Image("MaxwellVacuum/example_parfiles/Ex-convergence.png", width=500, height=500)


# In[26]:


Image("MaxwellVacuum/example_parfiles/Ax-convergence.png", width=500, height=500)


# <a id='errornodes'></a>
# 
# ## Step 5.b: Behavior of Error Nodes \[Back to [top](#toc)\]
# $$\label{errornodes}$$
# 
# Because System I is weakly hyperbolic (see [Tutorial-VacuumMaxwell_formulation_Cartesian](Tutorial-VacuumMaxwell_formulation_Cartesian.ipynb) for more discussion), zero speed error nodes of the constraint violation sit on our numerical grid, adding to the errors of our evolution variables. In contrast, System II is strongly hyperbolic, and the error nodes propagate away at the speed of light, leading to more stable evolution of the evolution variables. The plot below demostrates the qualitative behavior for both systems.
# 
# Contrast these plots to Figure 1 in [Knapp, Walker, & Baumgarte (2002)](https://arxiv.org/abs/gr-qc/0201051); we observe excellent qualitative agreement.

# In[27]:


Image("MaxwellVacuum/example_parfiles/constraintviolation.png", width=500, height=500)


# <a id='latex_pdf_output'></a>
# 
# # Step 6: Output this notebook to $\LaTeX$-formatted PDF file \[Back to [top](#toc)\]
# $$\label{latex_pdf_output}$$
# 
# The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename
# [Tutorial-ETK_thorn-MaxwellVacuum.pdf](Tutorial-ETK_thorn-MaxwellVacuum.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

# In[28]:


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