Implementation of the Green Operator/Discrete Green Operator hierarchy

[sbrisard]

Overview

Continuous Green operators

The (continuous) Green operator Γ is a linear operator that maps a field in ℝ^d onto a field in ℝ^d. It is attached to a specific physics: elasticity, conductivity, etc… In the remainder of this text, we will stick with the terminology of mechanics: the input of Γ will be referred to as the stress-polarization τ, while the output of Γ will be referred to as the strain ε. Note that the input and output can be tensors of any order; they are implemented in Scapin as column-vectors.

For periodic boundary conditions over the unit-cell Ω = (0, L₁) × (0, L₂) × …, the Green operator is defined in the Fourier space

Γ(τ)(x) =  Σ Γ^(kₙ)⋅τ^(kₙ) exp(ikₙ⋅x),
         n∈ℤ^d

where

          1  ⌠
τ^(k) =  ─── │ τ(x) exp(-ik⋅x) dx₁ dx₂ …
         |Ω| ⌡
             Ω

and

     ┌            ┐
     │ 2π n₁ / L₁ │
kₙ = │ 2π n₂ / L₂ │
     │ 2π n₃ / L₃ │
     └            ┘

The Fourier modes of Γ, Γ^(k), are known in closed form. They of course depend on the physics considered.

Discrete Green operators

For numerical simulations, the unit-cell Ω is discretized into a cartesian grid of size N = (N₁, N₂, …). The size of the cells is h = (h₁, h₂, …). 𝒫 denotes the set of cell-indices (in Julia language, 𝒫 = CartesianIndices(N)). The letters p and q will systematically be reserved to cell indices (in the real space), while the letters m and n will systematically be reserved to frequency indices (in the Fourier space).

We will use square brackets […] in the remainder for subscripting: τₕ[i, p₁, p₂, …] (or, in short, τₕ[i, p]) is the i-th component of the stress-polarization at cell p = (p₁, p₂, …). Likewise, τₕ^[i, n₁, n₂, …] (or, in short, τₕ^[i, n]) is the i-th component of the n-th mode of DFT(τₕ).

We need to consider discrete Green operators, hereafter denoted by a subscript “h” (Γₕ), which refers to the grid spacing. The discrete Green operator Γₕ is a linear operator that maps the discrete field τₕ (defined over the grid) onto a discrete field εₕ(also defined over the grid). All discrete Green operators implemented in Scapin are defined in the Fourier space

εₕ^[:, n] = Γₕ^[n]⋅τₕ^[:, n],

where τₕ^ and εₕ^ are the discrete Fourier transforms of τₕ and εₕ (along the last d axes). Therefore, application of the discrete Green operator reduces to the following operations: DFT of input → frequency-wise multiplication → inverse DFT.

The Fourier modes Γₕ^[n] of the discrete Green operator are generally defined as a linear combination of some Fourier modes of the continuous Green operator

Γₕ^[n] = Σ ϖ[𝓁] Γ^(q[n, 𝓁]),
         𝓁

where the scalars ϖ[𝓁] and vectors q[n, 𝓁] define the discretization scheme.

Implementation

In this section, we discuss various possible implementations of the Physical equations → Continuous Green operator → Discrete Green operator hierarchy of objects. Before we proceed, we observe that the following products

τ^ ↦ Γ^(k)⋅τ^    and    τₕ^ ↦ Γₕ^[n]⋅τₕ^

are implemented in a matrix-free form

apply(Γ^, k, τ^)    and    apply(Γₕ^, n, τₕ^),

the precise interface of the above functions is to be discussed in the present section.

Continuous Green operators

Option 1

Continuous Green operators are associated to a specific physical problem that must first be defined. For example

# Isotropic thermal conduction
struct Fourier{d,T}
    k::T
end

# Isotropic, linear elasticity
struct Hooke{d,T}
    μ::T
    ν::T
end

(d is the number of spatial dimensions, T is the scalar type). Then we can define the associated Green operators

struct GreenOperator{M}
    mat::M
end

where M is the material type (for lack of a better name) : Fourier, Hooke, …

Finally, we can define how Green operators operate in Fourier space

apply_fourier!(ε::AbstractVector{T}, Γ::GreenOperator{M}, τ::AbstractVector{T}, k::AbstractVector{T})
apply_fourier(Γ::GreenOperator{M}, τ::AbstractVector{T}, k::AbstractVector{T})

(equivalent of muladd should also be implemented). Note that GreenOperator should also have d and T in its list of type parameters.

Question 1: is apply_fourier a name that is precise enough? I'd rather use mul, which seems to be the Julia convention, but then I would have a name clash for discrete Green operators, where apply can operate both in the Fourier space and in the real space.

Question 2: should the signature be apply_fourier(Γ, τ, k) or apply_fourier(Γ, k, τ)? What does feel more natural? How would these methods be extended to allow for muladd-like operations?

Option 2

Objects of type Fourier and Hooke play a very minor role since GreenOperator{M} are really wrappers around M.

The benefit of the distinction between the material type and its associated Green operator is the fact that it is possible to define a apply operator for an object of type M: this object would apply the constitutive law to its input: σ = C * ε, where C is of type M.

The drawback is the multiplication of type parameters.

An alternative would be to directly define objects of type FourierGreenOperator, HookeGreenOperator, etc… and the apply_fourier!/apply_fourier methods.

Question 3: what is your opinion on option 2? Please elaborate!

Discrete Green operators

Discrete Green operators are associated to continuous Green operators:

struct DiscreteGreenOperator{G}
    Γ::G               # Underlying ContinuousGreenoperator
    N::Ntuple(d, Int)  # Grid size
    h::Ntuple(d, T)    # Grid spacing
end

Note that again, the underlying number of spatial dimensions d and scalar type T should be in the list of type parameters.

The τ^ ↦ Γₕ^[n]⋅τ^ operation would then be implemented as follows

apply_fourier!(εₕ::AbstractVector{T}, Γₕ::DiscreteGreenOperator{G}, τₕ::AbstractVector{T}, n::CartesianIndex{d})
apply_fourier(Γₕ::DiscreteGreenOperator{G}, τₕ::AbstractVector{T}, n::CartesianIndex{d})

And the τ^ ↦ Γ^(k)⋅τ^ convolution operation

apply!(εₕ::AbstractArray{T, d+1}, Γₕ::DiscreteGreenOperator{G}, τₕ::AbstractArray{T, d+1})
apply(Γₕ::DiscreteGreenOperator{G}, τₕ::AbstractArray{T, d+1})

which would be implemented in principle as follows (implementation is not syntactically correct and unoptimized in terms of memory allocation)

function apply(Γₕ::DiscreteGreenOperator{G}, τₕ::AbstractArray{T, d+1})
    τₕ^ = fft(τₕ, 2:(d+1))
    for n in CartesianIndices(Γₕ.N)
        apply_fourier!(view(εₕ^, :, n), Γ̂, n, τₕ^[:, n])
    end
    return ifft(εₕ^, 2:(d+1))
end

Question 4: is apply an appropriate function name? Should we prefer apply_real? Would it not be confusing, since sometimes, working in the real space would still require complex numbers (owing to time Fourier transforms – e.g. for viscoelasticity).

[antoine-levitt]

As a general scientific development principle, I tend to not think too much about APIs in advance. I know this is counter to what people say, that you should plan software development carefully, but I have just not found that to be useful, especially in scientific code whose purpose strongly depends on the research objectives of the moment. Indeed careful APIs easily lead to ossification of structure and then bloat. Rather, start from having something that works, even though it hurts aesthetics, and then play with it and figure out its pain points and the ways you might want to extend it later. We had several discussions about this with Michael at the beginning. He tends to be on the careful/planning side, and I tend to be on the careless/experimental one: we met in the middle and I think it was enriching for both of us (at least it was for me :-p)

That said I think your structure looks good? In general I would not sweat too much the structures: functions are more important than structures. Nothing forces you to use multiple dispatch if it's not natural for you.

If you wanted you could also use LinearOperators all the way, eg defining a RealToFourier and FourierToReal operator, and then doing something like (op::GreenOperator, x) = FourierToReal op.kernel RealToFourier x, with *(kernel, x) defined according to the physics. We have not done that in DFTK because we found functions to just be more flexible.

Regarding names we use real/fourier also even if sometimes the real quantities are complex numbers.

[sbrisard]

Thanks Antoine, I usually tend to aggre with you: start to write the program, then refactor. In the present case, I didn't even know where to start... but writing this post actually gave me a starting point!

That said I think your structure looks good? In general I would not sweat too much the structures: functions are more important than structures. Nothing forces you to use multiple dispatch if it's not natural for you.

I still have a lot to learn about Julia, because I don't really understand :'(

I am not sure about the general use of LinearOperators. I would not be able to use overloaded operators anyway, because I would like to operate in-place. Or did I miss something?

OK for real/fourier then.

Again, thanks for your time!

[antoine-levitt]

I still have a lot to learn about Julia, because I don't really understand :'(

I just meant that you don't necessarily have to fit in with julia's preferred way of doing things: the julia gods won't be mad at you if you don't use multiple dispatch and stick to a "C-like" imperative programming style. I did find that julia was not really friendly to OOP-style programming style, though.

I am not sure about the general use of LinearOperators. I would not be able to use overloaded operators anyway, because I would like to operate in-place. Or did I miss something?

I believe LinearMaps has support for in-place operations