Skip to content

Commit

Permalink
Merge pull request #120 from SciML/format
Browse files Browse the repository at this point in the history
Kind of format?
  • Loading branch information
ChrisRackauckas authored Feb 20, 2024
2 parents 8694b10 + ada91c3 commit 964deb0
Show file tree
Hide file tree
Showing 90 changed files with 2,347 additions and 1,866 deletions.
26 changes: 14 additions & 12 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
[![codecov](https://codecov.io/gh/SciML/MethodOfLines.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/SciML/MethodOfLines.jl)
[![Build Status](https://github.com/SciML/MethodOfLines.jl/workflows/CI/badge.svg)](https://github.com/SciML/MethodOfLines.jl/actions?query=workflow%3ACI)

[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
[![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)

MethodOfLines.jl is a package for automated finite difference discretization
Expand All @@ -19,16 +19,18 @@ This project is under active development, therefore the interface is subject to
Note that this package does not currently scale well to high resolution (high point count) problems, though there are changes in the works to remedy this.

Allowable terms in the system include, but are not limited to
- Advection
- Diffusion
- Reaction
- Nonlinear Diffusion
- Spherical Laplacian
- Any Julia function of the symbolic parameters/dependant variables and other parameters in the environment that's defined on the whole domain.

- Advection
- Diffusion
- Reaction
- Nonlinear Diffusion
- Spherical Laplacian
- Any Julia function of the symbolic parameters/dependant variables and other parameters in the environment that's defined on the whole domain.

Boundary conditions include, but are not limited to:
- Dirichlet
- Neumann (can also include time derivative)
- Robin (can also include time derivative)
- Periodic
- Any equation, which can include arbitrary Julia functions defined on that boundary, with the only symbolic parameters being those appearing in the referenced boundary.

- Dirichlet
- Neumann (can also include time derivative)
- Robin (can also include time derivative)
- Periodic
- Any equation, which can include arbitrary Julia functions defined on that boundary, with the only symbolic parameters being those appearing in the referenced boundary.
15 changes: 7 additions & 8 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,16 +6,15 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
# Make sure that plots don't throw a bunch of warnings / errors!
ENV["GKSwstype"] = "100"


include("pages.jl")

makedocs(sitename = "MethodOfLines.jl",
authors = "Chris Rackauckas, Alex Jones et al.",
clean = true, doctest = false, linkcheck = true,
modules = [MethodOfLines],
warnonly = [:docs_block, :missing_docs, :cross_references],
format = Documenter.HTML(assets = ["assets/favicon.ico"],
canonical = "https://docs.sciml.ai/MethodOfLines/stable/"),
pages = pages)
authors = "Chris Rackauckas, Alex Jones et al.",
clean = true, doctest = false, linkcheck = true,
modules = [MethodOfLines],
warnonly = [:docs_block, :missing_docs, :cross_references],
format = Documenter.HTML(assets = ["assets/favicon.ico"],
canonical = "https://docs.sciml.ai/MethodOfLines/stable/"),
pages = pages)

deploydocs(repo = "github.com/SciML/MethodOfLines.jl"; push_preview = true)
97 changes: 62 additions & 35 deletions docs/src/boundary_conditions.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
What follows is a set of allowable boundary conditions, please note that this is not exhaustive - try your condition and see if it works, the handling is quite general. If it doesn't, please post an issue and we'll try to support it. Currently, boundary conditions have to be supplied at the edge of the domain, but there are plans to support conditions embedded in the domain.

## Definitions

```julia
using ModelingToolkit, MethodOfLines, DomainSets

Expand All @@ -20,39 +21,48 @@ x_max = y_max = 1.0
```

## Dirichlet

```julia
v(t, 0, y) ~ 1.0
```

### Time dependent

```julia
u(t, 0., y) ~ x_min*y+ 0.5t
u(t, 0.0, y) ~ x_min * y + 0.5t
```

### Julia function

```julia
v(t, x, y_max) ~ sin(x)
```

### User-defined function

```julia
alpha = 9

f(t,x,y) = x*y - t
f(t, x, y) = x * y - t

function g(x,y)
z = sin(x*y)+cos(y)
function g(x, y)
z = sin(x * y) + cos(y)
# Note that symbolic conditionals require the use of IfElse.ifelse, or registration
return IfElse.ifelse(z > 0, x, 1.0)
end

u(t,x,y_min) ~ f(t,x,y_min) + alpha/g(x,y_min)
u(t, x, y_min) ~ f(t, x, y_min) + alpha / g(x, y_min)
```

### Registered User Defined Function

```julia
alpha = 9

f(t,x,y) = x*y - t
f(t, x, y) = x * y - t

function g(x,y)
z = sin(x*y)+cos(y)
function g(x, y)
z = sin(x * y) + cos(y)
# This function must be registered as it contains a symbolic conditional
if z > 0
return x
Expand All @@ -63,77 +73,94 @@ end

@register g(x, y)

u(t,x,y_min) ~ f(t,x,y_min) + alpha/g(x,y_min)
u(t, x, y_min) ~ f(t, x, y_min) + alpha / g(x, y_min)
```

## Neumann/Robin

```julia
v(t, x_min, y) ~ 2. * Dx(v(t, x_min, y))
v(t, x_min, y) ~ 2.0 * Dx(v(t, x_min, y))
```

### Time dependent

```julia
u(t, x_min, y) ~ x_min*Dy(v(t,x_min,y)) + 0.5t
u(t, x_min, y) ~ x_min * Dy(v(t, x_min, y)) + 0.5t
```

### Higher order

```julia
v(t, x, 1.0) ~ sin(x) + Dyy(v(t, x, y_max))
```

### Time derivative

```julia
Dt(u(t, x_min, y)) ~ 0.2
```

### User-defined function

```julia
function f(u, v)
(u + Dyy(v) - Dy(u))/(1 + v)
(u + Dyy(v) - Dy(u)) / (1 + v)
end

Dyy(u(t, x, y_min)) ~ f(u(t, x, y_min), v(t, x, y_min)) + 1
```

### 0 lhs

```julia
0 ~ u(t, x, y_max) - Dy(v(t, x, y_max))
```

## Periodic

```julia
u(t, x_min, y) ~ u(t, x_max, y)

v(t, x, y_max) ~ u(t, x_max, y)
```

Please note that if you want to use a periodic condition on a dimension with WENO schemes, please use a periodic condition on all variables in that dimension.

## Interfaces

You may want to connect regions with differing dynamics together, to do this follow the following example, splitting the variable that spans these domains:

```julia
@parameters t x1 x2
@variables c1(..)
@variables c2(..)
Dt = Differential(t)
@parameters t x1 x2
@variables c1(..)
@variables c2(..)
Dt = Differential(t)

Dx1 = Differential(x1)
Dxx1 = Dx1^2
Dx1 = Differential(x1)
Dxx1 = Dx1^2

Dx2 = Differential(x2)
Dxx2 = Dx2^2
Dx2 = Differential(x2)
Dxx2 = Dx2^2

D1(c) = 1 + c / 10
D2(c) = 1 / 10 + c / 10
D1(c) = 1 + c / 10
D2(c) = 1 / 10 + c / 10

eqs = [Dt(c1(t, x1)) ~ Dx1(D1(c1(t, x1)) * Dx1(c1(t, x1))),
Dt(c2(t, x2)) ~ Dx2(D2(c2(t, x2)) * Dx2(c2(t, x2)))]
eqs = [Dt(c1(t, x1)) ~ Dx1(D1(c1(t, x1)) * Dx1(c1(t, x1))),
Dt(c2(t, x2)) ~ Dx2(D2(c2(t, x2)) * Dx2(c2(t, x2)))]

bcs = [c1(0, x1) ~ 1 + cospi(2 * x1),
c2(0, x2) ~ 1 + cospi(2 * x2),
Dx1(c1(t, 0)) ~ 0,
c1(t, 0.5) ~ c2(t, 0.5), # Relevant interface boundary condition
-D1(c1(t, 0.5)) * Dx1(c1(t, 0.5)) ~ -D2(c2(t, 0.5)) * Dx2(c2(t, 0.5)), # Higher order interface condition
Dx2(c2(t, 1)) ~ 0]
bcs = [c1(0, x1) ~ 1 + cospi(2 * x1),
c2(0, x2) ~ 1 + cospi(2 * x2),
Dx1(c1(t, 0)) ~ 0,
c1(t, 0.5) ~ c2(t, 0.5), # Relevant interface boundary condition
-D1(c1(t, 0.5)) * Dx1(c1(t, 0.5)) ~ -D2(c2(t, 0.5)) * Dx2(c2(t, 0.5)), # Higher order interface condition
Dx2(c2(t, 1)) ~ 0]

domains = [t Interval(0.0, 0.15),
x1 Interval(0.0, 0.5),
x2 Interval(0.5, 1.0)]
domains = [t Interval(0.0, 0.15),
x1 Interval(0.0, 0.5),
x2 Interval(0.5, 1.0)]

@named pdesys = PDESystem(eqs, bcs, domains,
[t, x1, x2], [c1(t, x1), c2(t, x2)])
@named pdesys = PDESystem(eqs, bcs, domains,
[t, x1, x2], [c1(t, x1), c2(t, x2)])
```

Note that if you want to use a higher order interface condition, this may not work if you have no simple condition of the form `c1(t, 0.5) ~ c2(t, 0.5)`.
21 changes: 14 additions & 7 deletions docs/src/devnotes.md
Original file line number Diff line number Diff line change
@@ -1,14 +1,18 @@
# Notes for developers

## Getting started

First, fork the repo and clone it locally.

Then, type in the REPL

```
julia>] dev /path/to/your/repo
julia>] activate MethodOfLines
```

## Overview

MethodOfLines.jl makes heavy use of [`Symbolics.jl`](https://symbolics.juliasymbolics.org/dev/) and [`SymbolicUtils.jl`](https://symbolicutils.juliasymbolics.org), especially the replacement rules from the latter.

Take a look at [`src/discretization/MOL_discretization.jl`](https://github.com/SciML/MethodOfLines.jl/blob/master/src/MOL_discretization.jl) to get a high level overview of how the discretization works. A more consise description can be found [here](@ref hiw). Feel free to post an issue if you would like help understanding anything, or want to know developer opinions on the best way to go about implementing something.
Expand All @@ -21,13 +25,14 @@ A replacement rule is generated for each term which has a more specific higher s

Take a look at [`src/discretization/generate_finite_difference_rules.jl`](https://github.com/SciML/MethodOfLines.jl/blob/243252a595ed2af549d98270bd3b8ca5e3f93d69/src/discretization/generate_finite_difference_rules.jl) to see where the replacement rules are generated. Implemented schemes can be found in `/src/discretization/schemes`. Have a look at some of the already implemented examples there; read about the [`@rule` macro](https://symbolicutils.juliasymbolics.org/rewrite/) from `SymbolicUtils.jl`, if you haven't already. Note that the order that the rules are applied is important; there may be schemes that are applied first that are special cases of more general rules, for example the sphrical laplacian is a special case of the nonlinear laplacian.

First terms are split, isolating particular cases. Then, rules are generated and applied.
First terms are split, isolating particular cases. Then, rules are generated and applied.

Identify a rule which will match your case, then write a function that will handle how to apply that scheme for each index in the interior, for each combination of independant and dependant variables.
Identify a rule which will match your case, then write a function that will handle how to apply that scheme for each index in the interior, for each combination of independant and dependant variables.

This should be a function of the current index `II::CartesianIndex`, an independent variable `x` which represents the direction of the derivative, and a dependent variable `u`, which is the variable of which the derivative will be taken. The discrete representation of `u` is found in `s.discvars[u]`, which is an array with the same number of spatial dimensions as `u`, each index a symbol representing the discretized `u` at that index. Using this, and cartesian index offsets from `II`, create a finite difference/volume symbolic expression for the approximation of the derivative form you are trying to discretize. This should be returned.

For example, the following is a simple rule and function that would discretize derivatives of each dependent variable `u`in each dependent variable `x` with the second order central difference approximation:

```julia
#TODO: Add handling for cases where II is close to the boundaries
#TODO: Handle periodic boundary conditions
Expand All @@ -37,18 +42,20 @@ function second_order_central_difference(II::CartesianIndex, s::DiscreteSpace, u
j = x2i(s, u, x)

# Get a CartesianIndex of unit length that points in the direction of `x` e.g. CartesianIndex((1, 0, 0))
I1 = unitindex(ndims(u, s), j)
I1 = unitindex(ndims(u, s), j)

discu = s.discvars[u]
expr = (discu[II + I1] - discu[II - I1])/s.dx[x]
expr = (discu[II + I1] - discu[II - I1]) / s.dx[x]

return expr
end

# Note that indexmap is used along with the function `Idx` to create an equivalent index for the discrete form of `u`,
# which may have a different number of dimensions to `II`
function generate_central_difference_rules(II::CartesianIndex, s::DiscreteSpace, terms::Vector{<:Term}, indexmap::Dict)
rules = [[@rule Differential(x)(u) => second_order_central_difference(Idx(II, s, u, indexmap), s, u, x) for x in ivs(u, s)] for u in depvars]
function generate_central_difference_rules(
II::CartesianIndex, s::DiscreteSpace, terms::Vector{<:Term}, indexmap::Dict)
rules = [[@rule Differential(x)(u) => second_order_central_difference(
Idx(II, s, u, indexmap), s, u, x) for x in ivs(u, s)] for u in depvars]

rules = reduce(vcat, rules)

Expand All @@ -69,6 +76,6 @@ Initially, don't worry if your scheme is only implemented for specific approxima

Finally, include your rules in the vector of rules to be used to replace terms in the PDE at this index, found [here](https://github.com/SciML/MethodOfLines.jl/blob/949d0fee5e97c4adc59057460b3708161f776e9b/src/discretization/generate_finite_difference_rules.jl#L271):


## Inspecting generated code

To get the generated code for your system, use `code = ODEFunctionExpr(prob)`, or `MethodOfLines.generate_code(pdesys, discretization, "my_generated_code_filename.jl")`, which will create a file called `my_generated_code_filename.jl` in `pwd()`. This can be useful to find errors in the discretization, but note that it is not recommended to use this code directly, calling `solve(prob, AppropriateSolver())` will handle this for you.
Loading

0 comments on commit 964deb0

Please sign in to comment.