merging with main branch updates

.JuliaFormatter.toml

@@ -0,0 +1,2 @@
+style = "sciml"
+format_markdown = true

.github/workflows/CI.yml

@@ -3,9 +3,13 @@ on:
   pull_request:
     branches:
       - master
+    paths-ignore:
+      - 'docs/**'
   push:
     branches:
       - master
+    paths-ignore:
+      - 'docs/**'
 jobs:
   test:
     runs-on: ubuntu-latest
@@ -35,7 +39,7 @@ jobs:
           - '1'
           - '1.6'
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}

.github/workflows/Documentation.yml

@@ -11,7 +11,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'

.github/workflows/Downstream.yml

@@ -19,14 +19,14 @@ jobs:
         package:
           - {user: SciML, repo: PDESystemLibrary.jl, group: MOL}
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.julia-version }}
           arch: x64
       - uses: julia-actions/julia-buildpkg@latest
       - name: Clone Downstream
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
         with:
           repository: ${{ matrix.package.user }}/${{ matrix.package.repo }}
           path: downstream

.github/workflows/FormatCheck.yml

@@ -0,0 +1,42 @@
+name: format-check
+
+on:
+  push:
+    branches:
+      - 'master'
+      - 'release-'
+    tags: '*'
+  pull_request:
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: [1]
+        julia-arch: [x86]
+        os: [ubuntu-latest]
+    steps:
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: ${{ matrix.julia-version }}
+
+      - uses: actions/checkout@v4
+      - name: Install JuliaFormatter and format
+        # This will use the latest version by default but you can set the version like so:
+        #
+        # julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))'
+        run: |
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))'
+          julia  -e 'using JuliaFormatter; format(".", verbose=true)'
+      - name: Format check
+        run: |
+          julia -e '
+          out = Cmd(`git diff --name-only`) |> read |> String
+          if out == ""
+              exit(0)
+          else
+              @error "Some files have not been formatted !!!"
+              write(stdout, out)
+              exit(1)
+          end'

.github/workflows/Invalidations.yml

@@ -19,12 +19,12 @@ jobs:
     - uses: julia-actions/setup-julia@v1
       with:
         version: '1'
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - uses: julia-actions/julia-buildpkg@v1
     - uses: julia-actions/julia-invalidations@v1
       id: invs_pr
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         ref: ${{ github.event.repository.default_branch }}
     - uses: julia-actions/julia-buildpkg@v1

.gitignore

@@ -16,3 +16,4 @@ dev/*
 # Build artifacts for creating documentation generated by the Documenter package
 docs/build/
 docs/site/
+.DS_Store

Project.toml

@@ -1,7 +1,7 @@
 name = "MethodOfLines"
 uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
-authors = ["Alex Jones, <alex.lemony@gmail.com>"]
-version = "0.9.7"
+authors = ["Alex Jones, <alex.jones@juliahub.com>"]
+version = "0.10.2"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -15,6 +15,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PDEBase = "a7812802-0625-4b9e-961c-d332478797e5"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
@@ -30,10 +31,11 @@ Interpolations = "0.14"
 Latexify = "0.15, 0.16"
 ModelingToolkit = "8"
 OrdinaryDiffEq = "6"
-PDEBase = "0.1.4"
+PDEBase = "0.1.6"
 PrecompileTools = "1"
+RuntimeGeneratedFunctions = "0.5.9"
 SafeTestsets = "0.0.1"
-SciMLBase = "1.90, 1.91"
+SciMLBase = "1.94, 2"
 StaticArrays = "1"
 SymbolicUtils = "1"
 Symbolics = "5"
@@ -42,9 +44,10 @@ julia = "1.6"
 
 [extras]
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "NonlinearSolve", "SafeTestsets", "StableRNGs"]
+test = ["Test", "NonlinearSolve", "SafeTestsets", "StableRNGs", "Pkg"]

docs/Project.toml

@@ -10,10 +10,10 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
 DifferentialEquations = "7.6"
-Documenter = "0.27"
+Documenter = "1"
 DomainSets = "0.5, 0.6"
-MethodOfLines = "0.8, 0.9"
+MethodOfLines = "0.9, 0.10"
 ModelingToolkit = "8"
-NonlinearSolve = "1"
+NonlinearSolve = "1, 2"
 OrdinaryDiffEq = "6"
 Plots = "1"

docs/make.jl

@@ -10,19 +10,10 @@ include("pages.jl")
 
 makedocs(sitename = "MethodOfLines.jl",
          authors = "Chris Rackauckas, Alex Jones et al.",
-         clean = true,
-         doctest = false,
-         strict = [
-             :doctest,
-             :linkcheck,
-             :parse_error,
-             :example_block,
-             # Other available options are
-             # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block,
-         ],
+         clean = true, doctest = false, linkcheck = true,
          modules = [MethodOfLines],
-         format = Documenter.HTML(analytics = "UA-90474609-3",
-                                  assets = ["assets/favicon.ico"],
+         warnonly = [:docs_block, :missing_docs, :cross_references],
+         format = Documenter.HTML(assets = ["assets/favicon.ico"],
                                   canonical = "https://docs.sciml.ai/MethodOfLines/stable/"),
          pages = pages)
 

docs/pages.jl

@@ -9,6 +9,7 @@ pages = [
         "tutorials/icbc_sampled.md",
         "tutorials/PIDE.md",
     ],
+    "Performance Tips" => "performance.md",
     "MOLFiniteDifference" => "MOLFiniteDifference.md",
     "Solution Interface - PDESolutions" => "solutions.md",
     "Grid and Solution Retrieval - Deprecated" => "get_grid.md",

docs/src/howitworks.md

@@ -6,7 +6,7 @@ See [here](https://github.com/SciML/MethodOfLines.jl/blob/master/src/MOL_discret
 
 Given your discretization and `PDESystem`, we take each independent variable defined on the space to be discretized and create a corresponding range. We then take each dependent variable and create an array of symbolic variables to represent it in its discretized form. This is stored in a [`DiscreteSpace`](https://github.com/SciML/MethodOfLines.jl/blob/master/src/discretization/discretize_vars.jl) object, a useful abstraction.
 
-We recognize boundary conditions, i.e. whether they are on the upper or lower ends of the domain, or periodic [here](https://github.com/SciML/MethodOfLines.jl/blob/master/src/system_parsing/bcs/parse_boundaries.jl), and use this information to construct the interior of the domain for each equation [here](https://github.com/SciML/MethodOfLines.jl/blob/master/src/system_parsing/interiormap.jl). Each PDE is matched to each dependent variable in this step by which variable is of highest order in each PDE, with precedence given to time derivatives. This dictates which boundary conditions reduce the size of the interior for which PDE. This is done to ensure that there will be the same number of equations as discrete variable states, so that the system of equations is balanced.
+We recognize boundary conditions, i.e. whether they are on the upper or lower ends of the domain, or periodic [here](https://github.com/SciML/PDEBase.jl/blob/master/src/parse_boundaries.jl), and use this information to construct the interior of the domain for each equation [here](https://github.com/SciML/MethodOfLines.jl/blob/master/src/system_parsing/interior_map.jl). Each PDE is matched to each dependent variable in this step by which variable is of highest order in each PDE, with precedence given to time derivatives. This dictates which boundary conditions reduce the size of the interior for which PDE. This is done to ensure that there will be the same number of equations as discrete variable states, so that the system of equations is balanced.
 
 Next, the boundary conditions are discretized, creating an equation for each point on the boundary in terms of the discretized variables, replacing any space derivatives in the direction of the boundary with their upwind finite difference expressions. [This](https://github.com/SciML/MethodOfLines.jl/blob/master/src/discretization/generate_bc_eqs.jl) is the place to look to see how this happens.
 

docs/src/index.md

@@ -109,26 +109,19 @@ Pkg.status(;mode = PKGMODE_MANIFEST) # hide
 ```@raw html
 </details>
 ```
-```@raw html
-You can also download the 
-<a href="
-```
 ```@eval
 using TOML
-version = TOML.parse(read("../../Project.toml",String))["version"]
-name = TOML.parse(read("../../Project.toml",String))["name"]
-link = "https://github.com/SciML/"*name*".jl/tree/gh-pages/v"*version*"/assets/Manifest.toml"
-```
-```@raw html
-">manifest</a> file and the
-<a href="
-```
-```@eval
-using TOML
-version = TOML.parse(read("../../Project.toml",String))["version"]
-name = TOML.parse(read("../../Project.toml",String))["name"]
-link = "https://github.com/SciML/"*name*".jl/tree/gh-pages/v"*version*"/assets/Project.toml"
-```
-```@raw html
-">project</a> file.
+using Markdown
+version = TOML.parse(read("../../Project.toml", String))["version"]
+name = TOML.parse(read("../../Project.toml", String))["name"]
+link_manifest = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
+                "/assets/Manifest.toml"
+link_project = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
+               "/assets/Project.toml"
+Markdown.parse("""You can also download the
+[manifest]($link_manifest)
+file and the
+[project]($link_project)
+file.
+""")
 ```

docs/src/performance.md

@@ -0,0 +1,10 @@
+# Performance Tips
+
+While the default ModelingTolkit backend is highly capable, it doesn't scale well to very large systems of ODEs, present iun high resolution and/or higher dimensional PDE semidiscretizations that MethodOfLines.jl is capable of generating.
+
+To overcome this limitation, and scale MOL to realistic physics simulations, JuliaHub Inc. has released [JuliaSimCompiler.jl (installation instructions here)](https://help.juliahub.com/juliasimcompiler/dev/).
+
+Whenever this package and MethodOfLines are loaded together, MOL will also load `MethodOfLinesJuliaSimCompilerExt.jl`, an extension that takes advantage of this backend to speed up your compilation and solves automatically, with the same interface you are used to. 
+
+!!! note
+    JuliaSimCompiler is part of JuliaSim and thus requires a valid JuliaSim license to use. JuliaSim is a proprietary software developed by JuliaHub Inc. Using the packages through the registry requires a valid JuliaSim license. It is free to use for non-commercial academic teaching and research purposes. For commercial users, license fees apply. Please refer to the [End User License Agreement](https://juliahub.com/company/eula/?_gl=1*120lqg6*_ga*MTAxODQ4OTE3Mi4xNjk0MDA4MDM5*_ga_8FC7JQQLXX*MTY5NDAwODAzOC4xLjEuMTY5NDAwODgxMC4wLjAuMA..) for details. Please contact [[email protected]](https://juliahub.com/products/pricing/) for purchasing information.

docs/src/tutorials/heat.md

@@ -48,11 +48,11 @@ solu = sol[u(t, x)]
 using Plots
 plt = plot()
 
-for i in 1:length(discrete_t)
+for i in eachindex(discrete_t)
     plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])")
     scatter!(discrete_x, u_exact(discrete_x, discrete_t[i]), label="Exact, t=$(discrete_t[i])")
 end
-display(plt)
+plt
 ```
 ### Neumann boundary conditions
 
@@ -103,11 +103,11 @@ solu = sol[u(t, x)]
 using Plots
 plt = plot()
 
-for i in 1:length(discrete_t)
+for i in eachindex(discrete_t)
     plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])")
     scatter!(discrete_x, u_exact(discrete_x, discrete_t[i]), label="Exact, t=$(discrete_t[i])")
 end
-display(plt)
+plt
 ```
 
 ### Robin boundary conditions
@@ -159,11 +159,11 @@ solu = sol[u(t, x)]
 using Plots
 plt = plot()
 
-for i in 1:length(discrete_t)
+for i in eachindex(discrete_t)
     plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])")
     scatter!(discrete_x, u_exact(discrete_x, discrete_t[i]), label="Exact, t=$(discrete_t[i])")
 end
-display(plt)
+plt
 ```
 
 

src/MOL_discretization.jl

@@ -97,7 +97,7 @@ function SciMLBase.ODEFunction(pdesys::PDESystem, discretization::MethodOfLines.
                 gridlocs = get_gridloc.(us, (s,))
                 f_analytic = generate_function_from_gridlocs(analytic, gridlocs, s)
             end
-            return ODEFunction(simpsys; analytic = f_analytic, discretization.kwargs..., kwargs...)
+            return ODEFunction(simpsys; analytic = f_analytic, eval_module = @__MODULE__, discretization.kwargs..., kwargs...)
         end
     catch e
         println("The system of equations is:")

src/MethodOfLines.jl

@@ -14,6 +14,9 @@ using Interpolations
 using Latexify
 using PrecompileTools
 using DomainSets
+using RuntimeGeneratedFunctions
+RuntimeGeneratedFunctions.init(@__MODULE__)
+
 
 # See here for the main `symbolic_discretize` and `generate_system` functions
 using PDEBase
@@ -55,6 +58,7 @@ import SciMLBase.discretize
 # Interface
 include("interface/grid_types.jl")
 include("interface/scheme_types.jl")
+include("interface/callbacks.jl")
 include("interface/disc_strategy_types.jl")
 include("interface/MOLFiniteDifference.jl")
 
@@ -77,6 +81,7 @@ include("discretization/schemes/upwind_difference/upwind_diff_weights.jl")
 include("discretization/schemes/half_offset_weights.jl")
 include("discretization/schemes/extrapolation_weights.jl")
 include("discretization/differential_discretizer.jl")
+include("discretization/schemes/callbacks/callback_rules.jl")
 
 # System Parsing
 include("system_parsing/pde_system_transformation.jl")
@@ -110,6 +115,6 @@ include("precompile.jl")
 
 # Export
 export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete, chebyspace
-export UpwindScheme, WENOScheme, FunctionalScheme
+export UpwindScheme, WENOScheme, FunctionalScheme, MOLDiscCallback
 
 end

src/discretization/differential_discretizer.jl

@@ -8,13 +8,15 @@ struct DifferentialDiscretizer{T, D1, S} <: AbstractDifferentialDiscretizer
     interpmap::Dict{Num,DerivativeOperator}
     orders::Dict{Num,Vector{Int}}
     boundary::Dict{Num,DerivativeOperator}
+    callbacks
 end
 
 function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, discretization::MOLFiniteDifference, orders)
     pdeeqs = pdesys.eqs isa Vector ? pdesys.eqs : [pdesys.eqs]
     bcs = pdesys.bcs isa Vector ? pdesys.bcs : [pdesys.bcs]
     approx_order = discretization.approx_order
     advection_scheme = discretization.advection_scheme
+    callbacks = discretization.callbacks
     upwind_order = advection_scheme isa UpwindScheme ? advection_scheme.order : 1
 
     differentialmap = Array{Pair{Num,DerivativeOperator},1}()
@@ -63,5 +65,5 @@ function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, di
         push!(boundary, x => BoundaryInterpolatorExtrapolator(max(6, approx_order), dx))
     end
 
-    return DifferentialDiscretizer{eltype(orders),typeof(Dict(differentialmap)),typeof(advection_scheme)}(approx_order, advection_scheme, Dict(differentialmap), (Dict(nonlinlap_inner), Dict(nonlinlap_outer)), (Dict(windpos), Dict(windneg)), Dict(interp), Dict(orders), Dict(boundary))
+    return DifferentialDiscretizer{eltype(orders),typeof(Dict(differentialmap)),typeof(advection_scheme)}(approx_order, advection_scheme, Dict(differentialmap), (Dict(nonlinlap_inner), Dict(nonlinlap_outer)), (Dict(windpos), Dict(windneg)), Dict(interp), Dict(orders), Dict(boundary), callbacks)
 end

src/discretization/generate_finite_difference_rules.jl

@@ -80,8 +80,10 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace,
         integration_rules = []
     end
 
+    cb_rules = generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
+
     integration_rules = vcat(integration_rules, vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms)))
-    return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules)
+    return vcat(cb_rules, vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules)
 end
 
 function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace{W,M,G}, depvars, pde::Equation, derivweights::DifferentialDiscretizer, bmap, indexmap) where {W,M,G<:StaggeredGrid}

src/discretization/schemes/callbacks/callback_rules.jl

@@ -0,0 +1,5 @@
+function generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
+    cbs = derivweights.callbacks
+    cb_rules = [cb.sym => cb(s) for cb in cbs]
+    return cb_rules
+end