Numerical differential operators

CoorSystemNum provides the same differential operators as CoorSystemSym but evaluated pointwise at numerical coordinates using automatic differentiation (ForwardDiff.jl). This is useful when:

the field to differentiate is defined by a numerical function (not a symbolic expression);
a parametric study over many evaluation points is needed;
the coordinate system is defined only implicitly (e.g. from a position vector function).

Setup

Predefined factories return a CoorSystemNum ready to use:

using TensND, LinearAlgebra

CS_polar  = coorsys_polar_num()       # (r, θ)

CoorSystemNum{2, Float64}(TensND.var"#coorsys_polar_num##0#coorsys_polar_num##1"(), TensND.var"#coorsys_polar_num##2#coorsys_polar_num##3"(), TensND.var"#coorsys_polar_num##4#coorsys_polar_num##5"{TensND.var"#coorsys_polar_num##0#coorsys_polar_num##1"}(TensND.var"#coorsys_polar_num##0#coorsys_polar_num##1"()))

CS_cyl    = coorsys_cylindrical_num() # (r, θ, z)

CoorSystemNum{3, Float64}(TensND.var"#coorsys_cylindrical_num##0#coorsys_cylindrical_num##1"(), TensND.var"#coorsys_cylindrical_num##2#coorsys_cylindrical_num##3"(), TensND.var"#coorsys_cylindrical_num##4#coorsys_cylindrical_num##5"{TensND.var"#coorsys_cylindrical_num##0#coorsys_cylindrical_num##1"}(TensND.var"#coorsys_cylindrical_num##0#coorsys_cylindrical_num##1"()))

CS_sph    = coorsys_spherical_num()   # (θ, ϕ, r) — same convention as CoorSystemSym

CoorSystemNum{3, Float64}(TensND.var"#coorsys_spherical_num##0#coorsys_spherical_num##1"(), TensND.var"#coorsys_spherical_num##2#coorsys_spherical_num##3"(), TensND.var"#coorsys_spherical_num##4#coorsys_spherical_num##5"{TensND.var"#coorsys_spherical_num##0#coorsys_spherical_num##1"}(TensND.var"#coorsys_spherical_num##0#coorsys_spherical_num##1"()))

For a coordinate system defined by a position vector $\mathbf{OM}(x)$, use the generic constructor:

OM_func = x -> [x[1]*cos(x[2]), x[1]*sin(x[2]), x[3]];  # cylindrical
CS_gen  = CoorSystemNum(OM_func, 3)

CoorSystemNum{3, Float64}(TensND.var"#CoorSystemNum##0#CoorSystemNum##1"{Main.var"#2#3", Int64}(Main.var"#2#3"(), 3), TensND.var"#CoorSystemNum##4#CoorSystemNum##5"{Main.var"#2#3", Int64}(Main.var"#2#3"(), 3), TensND.var"#CoorSystemNum##8#CoorSystemNum##9"{Main.var"#2#3"}(Main.var"#2#3"()))

Operator syntax

All operators take a function and return a function. Evaluation at a point requires an extra call:

using TensND, LinearAlgebra
CS = coorsys_spherical_num();
f = x -> x[3]^2;      # f(θ,ϕ,r) = r²  in spherical coords
LAPLACE(f, CS)([π/4, π/3, 2.0])   # = n(n+1) r^(n-2) = 6 for n=2, r=2

6.0

The supported operators are:

Operator	Input	Output
`GRAD(f, CS)`	scalar → vector, vector → matrix	Function
`SYMGRAD(f, CS)`	vector → symmetric matrix	Function
`DIV(f, CS)`	vector → scalar, matrix → vector	Function
`LAPLACE(f, CS)`	scalar → scalar	Function
`HESS(f, CS)`	scalar → matrix	Function

All component arrays use the physical (normalized) basis matching the convention of CoorSystemSym.

Tensor results are returned as AbstractTens objects attached to the local normalized basis.

Basis accessors

The same accessors as CoorSystemSym are available, evaluated pointwise:

using TensND, LinearAlgebra
CS = coorsys_spherical_num();
x₀ = [π/3, π/4, 2.0];   # (θ, ϕ, r)
Lame(CS, x₀)             # Lamé coefficients χ = (r, r·sinθ, 1)

3-element Vector{Float64}:
 2.0
 1.7320508075688772
 1.0

normalized_basis(CS, x₀)

3×3 RotatedBasis{3, Float64}:
  0.353553  -0.707107  0.612372
  0.612372   0.707107  0.612372
 -0.866025   0.0       0.5

natvec(CS, x₀, 1, :cov)   # covariant natural vector aθ = χθ·eθ

3-element TensND.TensRotated{1, 3, Float64, Tensors.Vec{3, Float64}}:
 2.0
 0.0
 0.0

natvec(CS, x₀, 1, :cont)  # contravariant aθ = eθ/χθ

3-element TensND.TensRotated{1, 3, Float64, Tensors.Vec{3, Float64}}:
 0.5
 0.0
 0.0

unitvec(CS, x₀, 3)         # unit vector eʳ

3-element TensND.TensRotated{1, 3, Float64, Tensors.Vec{3, Float64}}:
 0.0
 0.0
 1.0

Basic scalar operators in polar coordinates

In polar coordinates $(r, \theta)$ with Lamé coefficients $(\chi_1, \chi_2) = (1, r)$, the classical result is

\[\nabla^2 (r^n f(\theta)) = (n^2 - m^2)\, r^{n-2} f(\theta) \qquad \text{for } f(\theta) = \sin(m\theta).\]

For $n = 2$, $m = 1$: $\nabla^2(r^2\sin\theta) = (4-1)\sin\theta = 3\sin\theta$.

using TensND
CS = coorsys_polar_num();
f = x -> x[1]^2 * sin(x[2]);   # r²sinθ
LAPLACE(f, CS)([2.5, 0.7])     # = 3·sin(0.7)

1.932653061713073

Gradient and Laplacian in spherical coordinates

For $f(r) = r^n$ in spherical coordinates $(\theta, \phi, r)$:

\[\nabla f = n r^{n-1}\, \mathbf{e}_r, \qquad \nabla^2 f = n(n+1)\, r^{n-2}.\]

For $f = 1/r$ (harmonic function): $\nabla f = -r^{-2}\mathbf{e}_r$, $\nabla^2 f = 0$.

using TensND, LinearAlgebra
CS = coorsys_spherical_num();
f = x -> x[3]^(-1);              # 1/r
GRAD(f, CS)([π/3, π/4, 2.0])    # → -1/r²·eʳ, components (eθ, eϕ, eʳ)

3-element TensND.TensRotated{1, 3, Float64, Tensors.Vec{3, Float64}}:
 -0.0
 -0.0
 -0.25

LAPLACE(f, CS)([π/3, π/4, 2.0]) # harmonic → 0

0.0

Lamé problem — hollow sphere under internal pressure

A hollow isotropic elastic sphere with inner radius $a$ and outer radius $b$ is subjected to an internal pressure $p$. The exact (Lamé) solution for the displacement is

\[u_r(r) = A\, r + \frac{B}{r^2}, \qquad A = \frac{p\,a^3}{3\kappa(b^3-a^3)}, \quad B = \frac{p\,a^3 b^3}{4\mu(b^3-a^3)},\]

and the non-zero strain components in the normalized spherical basis are

\[\varepsilon_{rr} = A - \frac{2B}{r^3}, \qquad \varepsilon_{\theta\theta} = \varepsilon_{\phi\phi} = A + \frac{B}{r^3}.\]

Parameters

using TensND, LinearAlgebra, Printf

CS = coorsys_spherical_num();   # coords: (θ, ϕ, r)
a, b   = 1.0, 3.0;              # inner / outer radius
p_i    = 1.0;                   # internal pressure
κ_val  = 2.0; μ_val = 1.0;     # bulk / shear modulus
λ_val  = κ_val - 2μ_val/3;
A = p_i * a^3 / (3κ_val * (b^3 - a^3))

0.00641025641025641

B = p_i * a^3 * b^3 / (4μ_val * (b^3 - a^3))

0.25961538461538464

Strain from SYMGRAD

The radial displacement field in the (θ, ϕ, r) normalized basis:

u_func = x -> [0.0, 0.0, A*x[3] + B/x[3]^2];
ε_func = x -> SYMGRAD(u_func, CS)(x);
ε = ε_func([π/3, π/4, 2.0])

3×3 TensND.TensRotated{2, 3, Float64, Tensors.SymmetricTensor{2, 3, Float64, 6}}:
 0.0388622  0.0         0.0
 0.0        0.0388622   0.0
 0.0        0.0        -0.0584936

The numerical values match the analytical expressions:

r₀ = 2.0;
A + B/r₀^3   # ε_θθ = ε_ϕϕ (exact)

0.03886217948717949

A - 2B/r₀^3  # ε_rr (exact)

-0.05849358974358975

Stress and equilibrium

The isotropic constitutive law $\boldsymbol{\sigma} = \lambda\operatorname{tr}(\boldsymbol{\varepsilon})\mathbf{1} + 2\mu\,\boldsymbol{\varepsilon}$:

σ_func = x -> begin
    ε    = Array(ε_func(x))
    tr_ε = sum(ε[i,i] for i in 1:3)
    [λ_val * tr_ε * (i==j ? 1.0 : 0.0) + 2μ_val * ε[i,j] for i in 1:3, j in 1:3]
end;

#10 (generic function with 1 method)

Verify equilibrium $\operatorname{div}\boldsymbol{\sigma} = \mathbf{0}$ and the boundary conditions $\sigma_{rr}(a) = -p$, $\sigma_{rr}(b) = 0$:

norm(Array(DIV(σ_func, CS)([π/3, π/4, 1.5])))  # ≈ 0 (machine precision)

2.2744152240584803e-17

σ_func([π/3, π/4, a])[3,3]   # σ_rr at inner radius ≈ -1

-1.0000000000000002

σ_func([π/3, π/4, b])[3,3]   # σ_rr at outer radius ≈ 0

-6.938893903907228e-18

Parametric study

println("  r     ε_rr(num)   ε_rr(exact)  σ_rr(num)  σ_rr(exact)  |div σ|")
for r in [1.1, 1.5, 2.0, 2.5, 2.9]
    x    = [π/3, π/4, r]
    ε    = Array(ε_func(x))
    σ    = σ_func(x)
    divσ = norm(Array(DIV(σ_func, CS)(x)))
    @printf("  %.1f  %+.6f  %+.6f   %+.6f  %+.6f   %.2e\n",
            r, ε[3,3], A-2B/r^3, σ[3,3], 3κ_val*A-4μ_val*B/r^3, divσ)
end

  r     ε_rr(num)   ε_rr(exact)  σ_rr(num)  σ_rr(exact)  |div σ|
┌ Warning: Assignment to `ε` in soft scope is ambiguous because a global variable by the same name exists: `ε` will be treated as a new local. Disambiguate by using `local ε` to suppress this warning or `global ε` to assign to the existing global variable.
└ @ coorsystems_num.md:216
  1.1  -0.383696  -0.383696   -0.741750  -0.741750   4.47e-16
  1.5  -0.147436  -0.147436   -0.269231  -0.269231   2.27e-17
  2.0  -0.058494  -0.058494   -0.091346  -0.091346   2.31e-18
  2.5  -0.026821  -0.026821   -0.028000  -0.028000   4.17e-17
  2.9  -0.014879  -0.014879   -0.004118  -0.004118   1.39e-17

Generic constructor from a position vector

For coordinate systems not predefined, pass an OM function to CoorSystemNum. The Lamé coefficients and Christoffel symbols are computed automatically by nested automatic differentiation of the metric tensor $g_{ij} = \partial_i \mathbf{OM} \cdot \partial_j \mathbf{OM}$.

Example — elliptic coordinates $(μ, ν)$ in 2D:

\[\mathbf{OM}(\mu, \nu) = a\begin{pmatrix}\cosh\mu\,\cos\nu \\ \sinh\mu\,\sin\nu\end{pmatrix}, \qquad \chi_1 = \chi_2 = a\sqrt{\sinh^2\!\mu + \sin^2\!\nu}.\]

using TensND
a_ell  = 2.0;
OM_ell = x -> [a_ell * cosh(x[1]) * cos(x[2]),
               a_ell * sinh(x[1]) * sin(x[2])];
CS_ell = CoorSystemNum(OM_ell, 2);
Lame(CS_ell, [1.0, 0.8])    # χ₁ = χ₂ = a√(sinh²μ + sin²ν)

2-element Vector{Float64}:
 2.75368669727873
 2.75368669727873

LAPLACE(x -> cosh(x[1])*cos(x[2]), CS_ell)([1.0, 0.8])  # harmonic → ≈ 0

1.387778780781445e-17