Forces

Different forces (or forcings, linear operators in general) are included in Limace.jl. For the high-level interface, they are implemented as abstract types Limace.Forcing{1} or Limace.Forcing{2}, for on one or two bases dependencies respectively.

Currently, the implemented forces are:

Some details for each of these forces is given in the following.

Inertial

\[\int \mathbf{u}_i^* \cdot \mathbf{u}_j\,\mathrm{d}V,\]

The linear operator is simply the inner product of the basis elements. For the momentum equation, this is essentially the projection of the inertial term $\partial_t \mathbf{u} \rightarrow \lambda \mathbf{u}$, hence the name.

Limace.Inertial — Type
mutable struct Inertial{TB, T} <: Limace.Forcing{1}
  • basis::Basis{TB}: The basis for the inertial operator.
  • factor::T: A scalar factor that multiplies the inertial operator, defaulting to 1.0.
  • mat::SparseMatrixCSC{ComplexF64}: A sparse matrix representation of the inertial operator.
  • preassembled::Bool: A flag indicating whether the inertial operator has been preassembled.

Example usage

u = Inviscid(10)
f = Limace.Inertial(u)
Limace.assemble!(f)
r.mat  # sparse matrix representation of the inertial operator
source
Limace.inertial — Function
inertial(
    b::Basis;
    threads,
    external
) -> LinearAlgebra.SymTridiagonal

Compute the Galerkin projection matrix of basis b onto itself, i.e. the inner products. Also known as the mass matrix.

source

Advection

\[\int \mathbf{u}_i^* \cdot \left(\left(\boldsymbol{\nabla}\times\mathbf{u}_j\right)\times\mathbf{U}_0 + \left(\boldsymbol{\nabla}\times\mathbf{U}_0\right)\times\mathbf{u}_j \right)\,\mathrm{d}V\]

This equivalent to the projection of the operator $\left(\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{U}_0+\left(\mathbf{U}_0\cdot\boldsymbol{\nabla}\right)\mathbf{u}$, when $\boldsymbol{\nabla}\cdot\mathbf{u} = \boldsymbol{\nabla}\cdot\mathbf{U}_0 = 0$.

Limace.Advection — Type
mutable struct Advection{TU, T} <: Limace.Forcing{1}
  • basis::Basis{TU}: The basis for the advection operator.
  • U0: The background flow, which can be a single BasisElement or a collection of them.
  • factor::T: A scalar factor that multiplies the advection operator, defaulting to 1.0.
  • mat::SparseMatrixCSC{ComplexF64}: A sparse matrix representation of the advection operator.
  • preassembled::Bool: A flag indicating whether the advection operator has been preassembled.

Example usage

u = Inviscid(10)
U0 = BasisElement(u, Toroidal, (1,0,0))
a = Limace.Advection(u, U0)
Limace.assemble!(a)
a.mat  # sparse matrix representation of the advection operator
source
Limace.advection — Function
advection(
    u::Basis,
    U0::BasisElement{T0<:Basis, TH<:Helmholtz, T};
    threads
) -> Any

Computes the advection term for a poloidal/toroidal background flow U0, a velocity basis u. Equivalent to -lorentz(u,u,U0).

source

Coriolis

\[\int \mathbf{u}_i^* \cdot 2\mathbf{e}_z\times\mathbf{u}_j\,\mathrm{d}V,\]

Limace.Coriolis — Type
mutable struct Coriolis{TB, T} <: Limace.Forcing{1}

  • basis::Basis

  • factor::Any

  • mat::SparseArrays.SparseMatrixCSC{ComplexF64}

  • preassembled::Bool

This defines the Coriolis force operator for a given basis. The factor is a scalar that multiplies the operator, e.g. the rotation rate $\Omega$ or a nondimensional parameter (e.g. $1/\mathrm{Le}$). The mat is a sparse matrix representation of the operator, and preassembled indicates whether the matrix has been preassembled.

Example usage

u = Inviscid(10)
Ω = 1.0
c = Limace.Coriolis(u, Ω)
Limace.assemble!(c)
c.mat  # sparse matrix representation of the Coriolis operator
source
Limace.coriolis — Function
coriolis(
    b::Basis;
    threads,
    Ω
) -> SparseArrays.SparseMatrixCSC{Tv, Int64} where Tv

Compute the sparse Galerkin projection matrix, by projecting the basis b onto the Coriolis operator.

source

Diffusion

\[\int \mathbf{u}_i^* \cdot \boldsymbol{\nabla}^2\mathbf{u}_j\,\mathrm{d}V\]

Limace.Diffusion — Type
mutable struct Diffusion{TB, T} <: Limace.Forcing{1}
  • basis::Basis{TB}: The basis for the diffusion operator.
  • factor::T: : A scalar factor that multiplies the diffusion operator, defaulting to 1.0.
  • mat::SparseMatrixCSC{ComplexF64}: A sparse matrix representation of the diffusion operator.
  • preassembled::Bool: A flag indicating whether the diffusion operator has been preassembled.
source
Limace.diffusion — Function
diffusion(
    b::Basis;
    threads,
    external
) -> SparseArrays.SparseMatrixCSC{Float64, Int64}

Compute the Galerkin projection matrix of the basis b onto the vector Laplacian. When keyword external=true, the integral is computed over all space, assuming continuity of the poloidal field and a scalar potential in the exterior domain.

source

Induction

\[\int \mathbf{b}_i^* \cdot \boldsymbol{\nabla}\times\left(\mathbf{u}_j\times\mathbf{B}_0\right)\,\mathrm{d}V\]

Limace.InductionB0 — Type
mutable struct InductionB0{TB, TU, T} <: Limace.Forcing{2}
  • bbasis::Basis{TB}: The magnetic field basis for the induction operator.
  • ubasis::Basis{TB}: The velocity basis for the induction operator.
  • B0: The background magnetic field, which can be a single BasisElement or a collection of them.
  • factor::T: : A scalar factor that multiplies the induction operator, defaulting to 1.0.
  • mat::SparseMatrixCSC{ComplexF64}: A sparse matrix representation of the induction operator.
  • preassembled::Bool: A flag indicating whether the induction operator has been preassembled.
u = Inviscid(10)
b = Insulating(10)
B0 = BasisElement(b, Poloidal, (2,0,1))
f = Limace.InductionB0(b, u, B0)
Limace.assemble!(f)
f.mat # sparse matrix representation of the induction operator
source

\[\int \mathbf{b}_i^* \cdot \boldsymbol{\nabla}\times\left(\mathbf{U}_0\times\mathbf{b}_j\right)\,\mathrm{d}V\]

Limace.InductionU0 — Type
mutable struct InductionU0{TB, T} <: Limace.Forcing{1}
  • basis::Basis{TB}: The basis for the induction operator.
  • U0: The background flow, which can be a single BasisElement or a collection of them.
  • factor::T: : A scalar factor that multiplies the induction operator, defaulting to 1.0.
  • mat::SparseMatrixCSC{ComplexF64}: A sparse matrix representation of the induction operator.
  • preassembled::Bool: A flag indicating whether the induction operator has been preassembled.

Example usage

u = Inviscid(10)
b = Insulating(10)
U0 = BasisElement(u, Toroidal, (2,0,1))
f = Limace.InductionU0(b, U0)
Limace.assemble!(f)
f.mat # sparse matrix representation of the induction operator
source
Limace.induction — Function
induction(
    bbi::Basis,
    buj::Basis,
    B0::BasisElement{T0<:Basis, TH<:Helmholtz, T};
    threads,
    external
) -> Any

Computes the induction term for a poloidal/toroidal background magnetic field B0, a magnetic field basis bbi and a velocity basis buj.

source
induction(
    bbi::Basis,
    U0::BasisElement{T0<:Basis, TH<:Helmholtz, T},
    bbj::Basis;
    threads,
    external
) -> Any

Computes the induction term for a poloidal/toroidal background velocity U0, a magnetic field basis bbi and a magnetic field basis bbj.

source

Lorentz

\[\int \mathbf{u}_i^* \cdot \left(\left(\boldsymbol{\nabla}\times\mathbf{b}_j\right)\times\mathbf{B}_0+\left(\boldsymbol{\nabla}\times\mathbf{B}_0\right)\times\mathbf{b}_j\right)\,\mathrm{d}V\]

Limace.Lorentz — Type
mutable struct Lorentz{TU, TB, T} <: Limace.Forcing{2}
  • ubasis::Basis{TU}: The velocity basis for the Lorentz operator.
  • bbasis::Basis{TB}: The basis for the Lorentz operator.
  • B0: The background magnetic field, which can be a single BasisElement or a collection of them.
  • factor::T: : A scalar factor that multiplies the Lorentz operator, defaulting to 1.0.
  • mat::SparseMatrixCSC{ComplexF64}: A sparse matrix representation of the Lorentz operator.
  • preassembled::Bool: A flag indicating whether the Lorentz operator has been preassembled.

Example usage

u = Inviscid(10)
b = Insulating(10)
B0 = BasisElement(b, Toroidal, (1,0,1))
f = Limace.Lorentz(u,b,B0)
Limace.assemble!(f)
f.mat = # sparse matrix representation of the Lorentz operator
source
Limace.lorentz — Function
lorentz(
    bui::Basis,
    bbj::Basis,
    B0::BasisElement{T0<:Basis, TH<:Helmholtz, T};
    threads
) -> Any

Computes the Lorentz term for a poloidal/toroidal background magnetic field B0, a velocity basis bui and a magnetic field basis bbj.

source

Explicit boundary condition operator

Limace.boundarycondition! — Method
boundarycondition!(A, b::TB, ::Type{T}=Float64) where {TB<:Basis,T<:Number}

Add the explicit boundary condition lines to a preassembled matrix A.

source
Limace.boundarycondition — Method
boundarycondition(b::TB, ::Type{T}=Float64) where {TB<:Basis,T<:Number}

Construct a sparse matrix representing the boundary conditions for the basis b.

source
Limace.zero_boundarycondition! — Method
zero_boundarycondition!(A, b::TB, ::Type{T}=Float64) where {TB<:Basis,T<:Number}

Zero out lines in a preassembled matrix (useful, when preassembling matrices at large resolution and then using them at lower resolution).

source

Assembly

All forcings are constructed without being assembled (i.e. the projections of each forcing onto the relevant basis are not done at the definition of the forcing).

The projections are done explicitly through Limace.assemble!, that saves the assembled sparse matrix in the forcing structure (f).

N = 10
u = Inviscid(N)
b = Insulating(N)
B0 = BasisElement(b, Poloidal, (1,0,1))
f = Limace.Lorentz(u,b,B0)
Limace.assemble!(f)
495×375 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2051 stored entries:
⎡⣀⠘⠦⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠣⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠉⢧⠀⠈⠓⢤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠙⣀⠀⠀⠹⢤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢦⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠓⣆⠀⠀⠙⠦⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠳⡄⠀⠀⠙⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢳⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠙⠢⣄⠀⠈⠙⢦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠹⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢳⡀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠲⢤⡀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠳⣄⎥
⎢⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠦⠀⠀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠳⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢤⡈⠓⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠘⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⡀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢩⡀⠀⠈⢲⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⡆⠀⠀⢳⣀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠣⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠘⢦⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢧⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡄⠀⠈⢳⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠳⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⣆⠀⠈⢇⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢷⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢳⡄⠈⢧⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠹⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢆⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠳⣄⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎦

This can also be done in parallel, by using the keyword argument threads=true.

Limace.assemble!(f; threads=true)

When the forcings are wrapped into a LimaceProblem, it is not necessary to call Limace.assembly! on each forcing by hand, but one can simply call Limace.assembly!(problem), as outlined in the Quickstart.