CalculustCore.Domains.BoxedDomain — Type
struct ProductDomain{T, D, var"#s10305"<:Tuple{Vararg{IntervalDomain, D}}} <: CalculustCore.Domains.AbstractDomain{T, D}

Logically rectangular product domain

CalculustCore.Domains.NullDomain — Type
struct NullDomain <: CalculustCore.Domains.AbstractDomain{Bool, -1}

Domain representing the empty set, ∅. Dimension set to -1.

CalculustCore.Domains.BoxDomain — Method
BoxDomain(endpoints; periodic_dims, tag, boundary_tags)

Creates a D-dimensional box domain with endpoints given by endpoints

Keyword Arguments

  • periodic_dims: List of periodic dimensions in 1:D. Defaults to none.
  • tag: Domain tag (symbol) defaults to :Interior.
  • boundary_tags: List of the 2D boundary tags. By default, periodic boundaries are tagged :Dk_periodic where k refers to the kth dimension, and aperiodic boundaries are tagged as :Dk_inf, :Dk_sup corresponding to the infimum and supremum in the kth dimension.

Example

julia> domain = BoxDomain(0, 1, 2, 3; periodic_dims=2)
(0, 1) × (2, 3)

julia> boundary_tags(domain)
(:D1_inf, :D1_sup, :D2_periodic, :D2_periodic)
CalculustCore.Domains.ChebyshevDomain — Method
ChebyshevDomain(D; kwargs...)

Creates a D-dimensional box domain with endpoints at -1.0, 1.0 in each dimension.

Keyword Arguments

  • periodic_dims: List of periodic dimensions in 1:D. Defaults to none.
  • tag: Domain tag (symbol) defaults to :Interior.
  • boundary_tags: List of the 2D boundary tags. By default, periodic boundaries are tagged :Dk_periodic where k refers to the kth dimension, and aperiodic boundaries are tagged as :Dk_inf, :Dk_sup corresponding to the infimum and supremum in the kth dimension.
CalculustCore.Domains.FixedEndpointBoxDomain — Method
FixedEndpointBoxDomain(D, x0, x1; kwargs...)

Creates a D-dimensional box domain with the same end points in each dimension.

Keyword Arguments

  • periodic_dims: List of periodic dimensions in 1:D. Defaults to none.
  • tag: Domain tag (symbol) defaults to :Interior.
  • boundary_tags: List of the 2D boundary tags. By default, periodic boundaries are tagged :Dk_periodic where k refers to the kth dimension, and aperiodic boundaries are tagged as :Dk_inf, :Dk_sup corresponding to the infimum and supremum in the kth dimension.
CalculustCore.Domains.FourierDomain — Method
FourierDomain(D; kwargs...)

Creates a D-dimensional box domain that is periodic in all directions. The domain spans from -π to π in each dimension.

Keyword Arguments

  • periodic_dims: List of periodic dimensions in 1:D. Defaults to none.
  • tag: Domain tag (symbol) defaults to :Interior.
  • boundary_tags: List of the 2D boundary tags. By default, periodic boundaries are tagged :Dk_periodic where k refers to the kth dimension, and aperiodic boundaries are tagged as :Dk_inf, :Dk_sup corresponding to the infimum and supremum in the kth dimension.
CalculustCore.Domains.UnitBoxDomain — Method
UnitBoxDomain(D; kwargs...)

Creates a D-dimensional, unit sized box domain with endpoints at 0.0, 1.0 in each dimension.

Keyword Arguments

  • periodic_dims: List of periodic dimensions in 1:D. Defaults to none.
  • tag: Domain tag (symbol) defaults to :Interior.
  • boundary_tags: List of the 2D boundary tags. By default, periodic boundaries are tagged :Dk_periodic where k refers to the kth dimension, and aperiodic boundaries are tagged as :Dk_inf, :Dk_sup corresponding to the infimum and supremum in the kth dimension.
CalculustCore.Domains.UnitCubeDomain — Method
UnitCubeDomain(; kwargs...)

Creates a unit-sized 3D cube domain spanning from 0.0 to 1.0 in each dimension.

Keyword Arguments

  • periodic_dims: List of periodic dimensions in 1:D. Defaults to none.
  • tag: Domain tag (symbol) defaults to :Interior.
  • boundary_tags: List of the 2D boundary tags. By default, periodic boundaries are tagged :Dk_periodic where k refers to the kth dimension, and aperiodic boundaries are tagged as :Dk_inf, :Dk_sup corresponding to the infimum and supremum in the kth dimension.
CalculustCore.Domains.UnitIntervalDomain — Method
UnitIntervalDomain(; kwargs...)

Creates a unit-sized 1D interval domain spanning from 0.0 to 1.0.

Keyword Arguments

  • periodic_dims: List of periodic dimensions in 1:D. Defaults to none.
  • tag: Domain tag (symbol) defaults to :Interior.
  • boundary_tags: List of the 2D boundary tags. By default, periodic boundaries are tagged :Dk_periodic where k refers to the kth dimension, and aperiodic boundaries are tagged as :Dk_inf, :Dk_sup corresponding to the infimum and supremum in the kth dimension.
CalculustCore.Domains.UnitSquareDomain — Method
UnitSquareDomain(; kwargs...)

Creates a unit-sized 2D square domain spanning from 0.0 to 1.0 in each dimension.

Keyword Arguments

  • periodic_dims: List of periodic dimensions in 1:D. Defaults to none.
  • tag: Domain tag (symbol) defaults to :Interior.
  • boundary_tags: List of the 2D boundary tags. By default, periodic boundaries are tagged :Dk_periodic where k refers to the kth dimension, and aperiodic boundaries are tagged as :Dk_inf, :Dk_sup corresponding to the infimum and supremum in the kth dimension.
CalculustCore.Domains.deform — Function
deform(domain; ...)
deform(domain, mapping; isseparable, isrescaling)

Deform domain via mapping, which expects the function signature

mapping(r1, ..., rD) = x1, ..., xD

where the inputs and outputs are NTuple{D, <:AbstractArray}.

Keyword Arguments

  • isrescaling - Set to true if mapping is a rescaling operation, i.e.

    x1 = a1 + λ1 × x1(r1),

    ...,

    xD = aD + λD × xD(rD)

  • isseparable - Set to true if mapping is separable, i.e.

    x1 = x1(r1),

    ...,

    xD = xD(rD)

TODO - make types for AffineMap, LinearMap, SeparableMap, Transation, Rotation TODO - check static_hasmethod in constructor

CalculustCore.Domains.expanse — Method
expanse(dom)

Returns size of domain bounding box.

Overload if expanse is not defined for bounding_box(dom).

CalculustCore.Spaces.Collocation — Type
struct Collocation <: CalculustCore.Spaces.AbstractDiscretization

A Collocation discretization is when the strong form of the equation is discretized as it is.

CalculustCore.Spaces.DeformedSpace — Type
struct DeformedSpace{T, D, Tspace<:CalculustCore.Spaces.AbstractSpace{T, D}, Tgrid, Tjacmat, Tjacimat, Tjac, Tjaci} <: CalculustCore.Spaces.AbstractSpace{T, D}

Deform domain, compute Jacobian of transformation, and its inverse

given 2D mapping

x = x(r,s), y = y(r,s)

compute

dXdR = [xr xs], dRdX = [rx ry], J = det(dXdR), Jinv = det(dRdX) [yr ys] [sx sy]

CalculustCore.Spaces.Galerkin — Type
struct Galerkin <: CalculustCore.Spaces.AbstractDiscretization

The Galerkin discretization is a weighted residual method where the trial space is equal to the test space.

Base.length — Method
length(V)

length of vector in space

Base.ndims — Method
ndims(_)

Dimension of underlying domain

Base.size — Method
size(V, d)

get number of points in dth dimension

CalculustCore.Spaces.advectionOp — Method
advectionOp(
    vels,
    V,
    discr;
    vel_update_funcs,
    vel_update_funcs!
)

Advection Operator: (v⃗⋅∇)⋅

for v,u,T in H¹₀(Ω)

(v,(u⃗⋅∇)T) = (v,ux∂xT + uy∂yT)

       = v' *B*(ux*∂xT + uy*∂yT)

implemented as

(u⃗⋅∇)T = ux∂xT + uy∂yT

   = [ux uy] * [Dx] T
               [Dx]
CalculustCore.Spaces.advectionOp — Method
advectionOp(vels, V, discr, Vd; ...)
advectionOp(vels, V, discr, Vd, J; vel_update_funcs)

When forming nonlinear operators...

CalculustCore.Spaces.advectionOp — Method
advectionOp(vels, V, discr, Vd; ...)
advectionOp(vels, V, discr, Vd, J; vel_update_funcs)

Galerkin discretizations suppport dealiased implementation of vector calculus operations where integration is performed in a(usually higher order) space, Vd. This is referred to as over-integration. If Vd is not provided, integration is done in V. J is the grid-to-grid interpolation operator from V to Vd. If J is not provided, it is recomputed. If Vd == V, then J becomes IdentityOperator(V).

ux,uy, ∇T are interpolated to a grid with higher polynomial order for dealiasing (over-integration) so we don't commit any "variational crimes"

CalculustCore.Spaces.biharmonicOp — Function
biharmonicOp(V::AbstractSpace, discr::AbstractDiscretization)

Biharmonic Operator: Δ²

Represents the Biharmonic opeator over V per discretization scheme discr.

CalculustCore.Spaces.diffusionOp — Method
diffusionOp(ν::Number, V::AbstractSpace, discr::AbstractDiscretization)

Diffusion operator: -νΔ

Returns the negative Laplace Operator scaled by diffusion coefficient ν. ν can be updated per ν_update_func which expects the same signature as SciMLOperators.ScalarOperator. All positional arguments besides ν are passed down to laplaceOp to form the negative lapalcian.

CalculustCore.Spaces.diffusionOp — Method
diffusionOp(ν, V, discr; ν_update_func)

Diffusion operator: -∇⋅(ν∇⋅)

Returns the diffusion operator where ν is the space-varying diffusion coefficient in V. ν can be updated per ν_update_func which expets the same signature as SciMLOperators.DiagonalOperator.

It is assumed that ν is a function in space V. However, if V is a transformed space, then ν is a function in the physical space, i.e. transform(V).

CalculustCore.Spaces.diffusionOp — Method
diffusionOp(ν, V, discr, Vd; ...)
diffusionOp(ν, V, discr, Vd, J; ν_update_func)

Galerkin discretizations suppport dealiased implementation of vector calculus operations where integration is performed in a(usually higher order) space, Vd. This is referred to as over-integration. If Vd is not provided, integration is done in V. J is the grid-to-grid interpolation operator from V to Vd. If J is not provided, it is recomputed. If Vd == V, then J becomes IdentityOperator(V).

CalculustCore.Spaces.forcingOp — Method
forcingOp(f, V, discr; f_update_func, f_update_func!)

Added forcing as an operator.

F = forcingOp(f) L = A + F

F(u) = 0 * u + M * f L(u) = A * u + M * f

CalculustCore.Spaces.form_transform — Function
form_transform(V::AbstractSpace, input::AbstractVecOrMat; [kwargs...])

Form transform operator given by transformOp(V) that may be applied to the provided input prototype array per keyword arguments kwargs. kwargs may include p, the parameter set.

CalculustCore.Spaces.gradientOp — Method
gradientOp(V::AbstractSpace{T, D}) where{T, D} -> [∂x1, ..., ∂xD]

Gradient Operator: ∇

Returns a size D Vector of operators. The dth operator corresponds to the linear transformation from a function to its partial derivative in the dth dimension at the grid points of V. Its [ij]th entry is equal to the value of the derivative of the jth basis function at the ith grid point.

$[D]_{ij} = \partial_{x}\phi_{j}(x_i)$

CalculustCore.Spaces.hessianOp — Method
hessianOp(V::AbstractSpace{T, D}) where{T, D} ->

$[∂x1∂x1 ... ∂x1∂xD, ... ∂xD∂x1 ... ∂xD∂xD ]$

Hessian Operator: ∇²

Returns a D × D Matrix of operators whose [ij]th entry represents the transformation from a function to its second partial derivative in the [ij]th directions.

CalculustCore.Spaces.interp — Function
interp(points, u::AbstractVector, V::AbstractSpace)

Interpolate function described by u to points points.

CalculustCore.Spaces.interpOp — Function
interpOp(V1::AbstractSpace, V2::AbstractSpace)

Grid to grid interpolation operator from V1 to V2. Both spaces must share the same domain. Its [ij]th entry is the value of the jth basis function of V1 at the ith grid point of V2.

$[J]_{ij} = \phi_{j}(x_i)$

If V1 == V2, we simply return the no-op IdentityOperator(V1).

CalculustCore.Spaces.laplaceOp — Function
laplaceOp(V::AbstractSpace, discr::AbstractDiscretization)

Laplace Operator: -Δ

Returns the negative Laplace operator over V per discretization scheme discr. As the Laplace operator (laplacian) is a negative-definitie operator, we return the negative laplacian which is a positive-definite operator.

CalculustCore.Spaces.laplaceOp — Method

(v,-∇² u) = (vx,ux) + (vy,uy)

     := a(v,u)

      = v' * A * u

      = (Q*R'*v)'*A_l*(Q*R'*u)

      = v'*R*Q'*A_l*Q*R'*u

implemented as

R'R * QQ' * Al * uloc where A_l is

[Dr]'[rx sx]'[B 0][rx sx][Dr] [Ds] [ry sy] [0 B] [ry sy] [Ds]

= [Dr]' * [G11 G12]' * [Dr] [Ds] [G12 G22] [Ds]

CalculustCore.Spaces.laplaceOp — Method
laplaceOp(V, discr, Vd)
laplaceOp(V, discr, Vd, J)

Galerkin discretizations suppport dealiased implementation of vector calculus operations where integration is performed in a(usually higher order) space, Vd. This is referred to as over-integration. If Vd is not provided, integration is done in V. J is the grid-to-grid interpolation operator from V to Vd. If J is not provided, it is recomputed. If Vd == V, then J becomes IdentityOperator(V).

CalculustCore.Spaces.laplaceOp — Method
laplaceOp(V, discr)

For a Galerkin discretization, the entries of the negative laplacian are given by

$[A]_{ij} = \int_\Omega \nabla\phi_i(x) \cdot \nabla\phi_j(x) dx$

for v,u in H¹₀(Ω)

(v,-∇² u) = (vx,ux) + (vy,uy)

     := a(v,u)
CalculustCore.Spaces.make_transform — Function
make_transform(V; ...)
make_transform(V, input; kwargs...)

Returns a copy of V with transform operator (given by transformOp(V)) that may be applied to the provided input prototype array per provided keyword arguments kwargs.

CalculustCore.Spaces.massOp — Function
massOp(V::AbstractSpace, discr::AbstractDiscretization)

Mass Operator: ∫⋅dx

Represents the mass matrix for discretization scheme discr.

For a Galerkin discretization, the [ij]th entry corresponds to the inner product of the ith and jth basis functions.

$[M]_{ij} = \int_\Omega \phi_i(x)\phi_j(x) dx = \langle \phi_i, \phi_j \rangle$

CalculustCore.Spaces.massOp — Method
massOp(V, discr, Vd)
massOp(V, discr, Vd, J)

Galerkin discretizations suppport dealiased implementation of vector calculus operations where integration is performed in a(usually higher order) space, Vd. This is referred to as over-integration. If Vd is not provided, integration is done in V. J is the grid-to-grid interpolation operator from V to Vd. If J is not provided, it is recomputed. If Vd == V, then J becomes IdentityOperator(V).

CalculustCore.Spaces.transformOp — Function
transformOp(V::AbstractSpace)

Get transform operator that takes vectors from V to transform(V). To change transform operator to act on a different AbstractVecOrMat subtype, call

V = make_transform(V, input_prototype)
CalculustCore.Spaces.truncationOp — Function
truncationOp(V::AbstractSpace{T, D} , fracs::NTuple{D}) where{T, D}

Truncation operator removes high-frequency modes in an input vector. fracs[d] is the proprtion of the spectrum to be preserved in dth dimension. In spectral space, the operator corresponds to diagonal scaling, where the first fracs[d] of the spectrum (low frequency modes) is multiplied by unity, and the remainder (high-frequency modes) with zeros.

CalculustCore.CalculustCore — Module

CalculustCore.jl

A common interface for PDE solving

PDEs are messy: weird boundary conditions, moving domains, singular operators, you name it! And solvers more so with their discretization schemes, stability charectoristics, and performance tradeoffs. Not to mention the layers of complications added by the software stack, or computer architecture, or pre/post-processing pipelines.

Numerical PDE problems drive high-performance computing. The biggest supercomputers on the planet run physics simulations on millions of MPI ranks over years of wall-time. An optimized workflow is the difference between having a solution, or going home emptyhanded, and every ounce of performance has to be squeezed. As such, highly tuned software packages that specialize on a set class of problems dominate the market. With specializations, however, generalizability and interpoeratibility take a hit.

CalculustCore.jl is written so that package authors won't need to write a method for Δ (Laplacian) every time a new scheme for the gradient ∇ (gradient) comes along. We provides abstractions over components of PDE solvers that reduce the amount of boilerplate code that any new solver would need, and improve interpoeratibility between solvers. Furthermore, a requisite for the developement of Machine Learning based PDE solvers is the ability to mix state of the art discretizations with large neural network models on a AD-compatible, accelerator-friendly test-bench. Julia's multiple-dispatch based programming model allows for abstractly-typed, composable code, that lets packages just plug-and-play. Engineers shouldn't waste their time reinventing the wheel, rather focus on what's new and interesting!

We believe switching from 'Fourier-spectral collocation method' to 'discontinuous Galerkin finite elements' should be a one-line change for the user. And users should not need to deal with inconsistent syntax between solvers for specifying boundary conditions, and forming vector-calculus operators.

Finally, CalculustCore.jl is fully compatible with the SciML.ai ecosystem, so after describing your problem, it should spit out the right BoundaryValuePDEProblem or ODEProblem that you can solve using the relevant DiffEq package.

Abstract Interfaces

CalculustCore.jl contains separate abstract interaces for multidimensional domains, vector-calculus operators, and function spaces. It is general enough that anybody can plug in their discretizations (ie define an inner-product operator, and gradient operator) and start solving boundary value problems, or time-evolution problems.

Once you plug in your discretizations, you can do a lot of cool things like apply any random deformations to the space. CalculustCore.jl translate all your vector calculus operators correctly. That means the same code could solve convection diffusion on a square as well as an annulus with no extra work and basically conserved accuracy.

AbstractDomain interface

concrete types, Boundary tags

  • deform - all the mapping stuff
  • \otimes

concrete types

  • IntervalDomain
  • BoxDomain

AbstractSpace interface

define these methods (grad, mass, common functions)

  • deform
  • \otimes
  • transform

AbstractDiscretization interface

  • GalerkinProjection
  • Collocation

Operator interface - SciMLOperators.jl

BoundaryValueProblem interface

Usually a rank-deficient systems

Associated Packages

Roadmap

  • [ ] AbstractDomain interface

    • [ ] Move concrete types to a separate package - PDEDomains.jl
    • [X] Logically rectangular domains
    • [X] Deform domain
    • [X] Boundary tags
    • [ ] Interior tags (for multiphase flows, conjugate heat-transfer)
    • [ ] Gordon Hall interpolation (transfinite interpolation)
    • [ ] meshed domains
    • [ ] signed distance geometries
    • [ ] Time-varying domains
    • [ ] Is it possibe to just use DomainSets.jl and add some metadata info?

  • [ ] AbstractField <: AbstractVector interface - special array types

    • [X] Spectral polynomial (nothing special needed)
    • [X] transform-based spectral (fourier, cheby) (nothing special needed)
    • [ ] Box/ full spectral elements

  • [ ] Operator interface - moved to SciMLOperators.jl

    • [X] linear algebra operations
    • [X] lazy composition
    • [X] can use array reductions
    • [X] caching
    • [ ] Gather-Scatter operator using NNlib
    • [ ] General interpolation operator on element-meshes

  • [ ] Spaces

    • [X] Deformed spaces

    • [ ] Tensor product spaces

    • [X] transformed space

    • [X] orthogonal polynomials - NodalPolynomialSpaces.jl

    • [ ] Spectral with transforms

      • [X] Fourier - FourierSpaces.jl
      • [ ] Cosine/ Sine spaces

    • [ ] Box, full spectral elements - create SpectralElementSpaces.jl

  • [X] Create a distinction between Space, and Discretization

    • [X] Space is how to represent functions
    • [X] Discretization is how you form operators
    • [ ] Rename "discretization" to "scheme" or "solve_scheme" or something. Because discretization is ambiguous. is it referring to spatial discretization? or time discretization? In this package we are using it as a "scheme" to form differential operators on a "discretized" space.
    • [ ] Rename Galerkin -> GalerkinProjection for clarity
    • [ ] Flux handling in discontinuous galerkin/ finite volume
    • [ ] all the shennanigans in stabalized finite elements

  • [ ] Boundary Condition interface (apply "this" boundary condition based on "that" domain boundary)

    • [X] Dirichlet
    • [X] Neumann
    • [ ] Robin

  • [ ] Problems

    • [ ] Problem frontend with ModelingToolkit.jl

    • [ ] Boundary Value Problems

      • [ ] move boundary information to RHS
      • [ ] dispatch to LinearSolve.jl
      • [ ] dispatch to NonlinearSolve.jl (after LinearSolve.jl, NonlinearSolve.jl integration)

    • [ ] Time-evolution problems

      • [X] play nice with OrdinaryDiffEq
      • [ ] for implicit time-steppers, solve a BVP at every time step. impose boundary condition on the operator (wait for SciMLOperators, OrdinaryDiffEq integration)
      • [ ] automatically spit out an ODEProblem