Set up ghost pixels (#12)

Project.toml

@@ -15,6 +15,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

examples/dipole.jl

@@ -23,13 +23,13 @@ sources = [
 sim = Khronos.Simulation(
     cell_size = [4.0, 4.0, 4.0],
     cell_center = [0.0, 0.0, 0.0],
-    resolution = 40,
+    resolution = 10,#40,
     sources = sources,
     boundaries = [[1.0, 1.0], [1.0, 1.0], [1.0, 1.0]],
 )
 
 # Run the simulation for 10 "time units"
-t_end = 40.0;
+t_end = 1.0#40.0;
 Khronos.run(sim, until = t_end)
 
 # Plot cross sections of the result

examples/planewave.jl

@@ -5,6 +5,7 @@
 import Khronos
 using CairoMakie
 using GeometryPrimitives
+using LinearAlgebra
 
 Khronos.choose_backend(Khronos.CUDADevice(), Float64)
 
@@ -23,13 +24,14 @@ function build_simulation(λ, resolution)
     cell_size = [sxy, sxy, sz]
 
     boundary_layers = [[dpml, dpml], [dpml, dpml], [dpml, dpml]]
-
+    k_vector = [0.00, 0.0, 1.0]
+    k_vector = k_vector ./ norm(k_vector)
     sources = [
         Khronos.PlaneWaveSource(
             time_profile = Khronos.ContinuousWaveSource(fcen = 1.0 / λ),
             center = [0, 0, 0],
             size = [Inf, Inf, 0],
-            k_vector = [0.0, 0.0, 1.0],
+            k_vector = k_vector,
             polarization_angle = 0.0,
             amplitude = 1im,
         ),

src/DataStructures.jl

@@ -326,7 +326,11 @@ function get_Δz(cell_size, Nz)
     return cell_size[3] / (Nz)
 end
 
-@with_kw mutable struct SimulationData{N,T,CN,CT}
+@with_kw mutable struct SimulationData{N,T,CN,CT,BT}
+    # N: Real number type
+    # T: Real array type
+    # CN: Complex number type
+    # CT: Complex array type
     cell_size::Vector{N} = N.(cell_size)
     cell_center::Any
     resolution::Any
@@ -348,7 +352,7 @@ end
     dimensionality::Type{<:Dimension} = ThreeD
 
     # Data to be generated upon initialization
-    fields::Union{Fields{T},Nothing} = nothing
+    fields::Union{Fields{BT},Nothing} = nothing
     source_data::Union{Vector{SourceData{CT}},Nothing} = nothing
     source_components::Union{Vector{<:Field},Nothing} = nothing
     boundary_data::Union{BoundaryData{T},Nothing} = nothing
@@ -363,6 +367,7 @@ Simulation(; kwargs...) = SimulationData{
     backend_array,
     complex_backend_number,
     complex_backend_array,
+    AbstractArray,
 }(;
     kwargs...,
 )

src/Fields.jl

@@ -1,5 +1,9 @@
 # (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary.
 
+# -------------------------------------------------------------------------- #
+# General field array utility functions
+# -------------------------------------------------------------------------- #
+
 """
     next_dir(direction)
 
@@ -89,93 +93,141 @@ get_fields_from_component(sim::SimulationData, ::Hz) = sim.fields.fHz
 # FIXME: once we have offdiag compoents, the inverse is more complicated...
 get_fields_from_component(sim::SimulationData, ::εx) = inv.(sim.geometry_data.ε_inv_x)
 get_fields_from_component(sim::SimulationData, ::εy) = inv.(sim.geometry_data.ε_inv_y)
-get_fields_from_component(sim::SimulationData, ::εz) = inv.(sim.geometry_data.ε_inv_z)
+function get_fields_from_component(sim::SimulationData, ::εz)
+    inv.(sim.geometry_data.ε_inv_z)
+end
+
+# -------------------------------------------------------------------------- #
+# General simulation field allocation functions.
+#
+#
+# [smartalecH] Currently, Khronos.jl only supports full 3D simulations. In the
+# future, however, it would be nice to support 2D simulations (and
+# cylindrically-symmetric simulations). In this case, one could specify sources
+# and materials that only require half as many computational resources as a full
+# 2D simulation. Specifically, they only require a TE or a TM configuration. We
+# can create independent functions that set up field arrays for these disjoint
+# cases, and then rely on dispatch to pick the proper configuration once Khronos
+# is ready for it. Since 3D is just a superset of the two configurations, we'll
+# just call both.
+# -------------------------------------------------------------------------- #
+
+"""
+    initialize_field_component_array(sim::SimulationData, component::Field)::AbstractArray
+
+Initialize an electromagnetic field array for a given `component` over the
+entire simulation (`sim`) domain.
+
+We add a "ghost" pixel to the boundary of our domain to more easily facilitate
+boundary conditions with the GPU. In the future, we can use this same formalism
+to enable MPI communication between "chunks" of data.
+
+To ensure that these ghost pixels serve as a "drop-in" with the existing stencils,
+we use offset arrays and start 1-based indexing after the first ghost pixel.
+"""
+function initialize_field_component_array(
+    sim::SimulationData,
+    component::Field,
+)::AbstractArray
+    array_dims = get_component_voxel_count(sim, component) .+ 2
+    return OffsetArray(
+        KernelAbstractions.zeros(backend_engine, backend_number, array_dims...),
+        -1,
+        -1,
+        -1,
+    )
+end
 
 """
-    get_fields_args_TE(sim)
+    allocate_fields_TE(sim)
 
-Pull the proper fields arguments for a 2D TE simulation.
+Allocate the proper field components for a 2D TE simulation (where the E-field is
+"out-of-plane"). This includes Dz/Ez, Bx/Hx, By/Hy, and their corresponding
+current sources. 
 """
-function get_fields_args_TE(sim::SimulationData)
+function allocate_fields_TE(sim::SimulationData)
     return (
-        fEz = zeros(get_component_voxel_count(sim, Ez())[1:sim.ndims]...),
-        fDz = zeros(get_component_voxel_count(sim, Dz())[1:sim.ndims]...),
-        fCDz = needs_C(sim, Dz()) ?
-               zeros(get_component_voxel_count(sim, Dz())[1:sim.ndims]...) : nothing,
-        fUDz = needs_U(sim, Dz()) ?
-               zeros(get_component_voxel_count(sim, Dz())[1:sim.ndims]...) : nothing,
-        fWDz = needs_W(sim, Dz()) ?
-               zeros(get_component_voxel_count(sim, Dz())[1:sim.ndims]...) : nothing,
-        fHx = zeros(get_component_voxel_count(sim, Hx())[1:sim.ndims]...),
-        fBx = zeros(get_component_voxel_count(sim, Bx())[1:sim.ndims]...),
-        fCBx = needs_C(sim, Bx()) ?
-               zeros(get_component_voxel_count(sim, Bx())[1:sim.ndims]...) : nothing,
-        fUBx = needs_U(sim, Bx()) ?
-               zeros(get_component_voxel_count(sim, Bx())[1:sim.ndims]...) : nothing,
-        fWBx = needs_W(sim, Bx()) ?
-               zeros(get_component_voxel_count(sim, Bx())[1:sim.ndims]...) : nothing,
-        fHy = zeros(get_component_voxel_count(sim, Hy())[1:sim.ndims]...),
-        fBy = zeros(get_component_voxel_count(sim, By())[1:sim.ndims]...),
-        fCBy = needs_C(sim, By()) ?
-               zeros(get_component_voxel_count(sim, By())[1:sim.ndims]...) : nothing,
-        fUBy = needs_U(sim, By()) ?
-               zeros(get_component_voxel_count(sim, By())[1:sim.ndims]...) : nothing,
-        fWBy = needs_W(sim, By()) ?
-               zeros(get_component_voxel_count(sim, By())[1:sim.ndims]...) : nothing,
-        fSBx = Hx() ∈ sim.source_components ?
-               zeros(get_component_voxel_count(sim, Hx())[1:sim.ndims]...) : nothing,
-        fSBy = Hy() ∈ sim.source_components ?
-               zeros(get_component_voxel_count(sim, Hy())[1:sim.ndims]...) : nothing,
-        fSDz = Ez() ∈ sim.source_components ?
-               zeros(get_component_voxel_count(sim, Ez())[1:sim.ndims]...) : nothing,
+        fEz = initialize_field_component_array(sim, Ez()),
+        fDz = initialize_field_component_array(sim, Dz()),
+        # Only allocate PML auxilliary fields if needed
+        fCDz = needs_C(sim, Dz()) ? initialize_field_component_array(sim, Dz()) : nothing,
+        fUDz = needs_U(sim, Dz()) ? initialize_field_component_array(sim, Dz()) : nothing,
+        fWDz = needs_W(sim, Dz()) ? initialize_field_component_array(sim, Dz()) : nothing,
+        fHx = initialize_field_component_array(sim, Hx()),
+        fBx = initialize_field_component_array(sim, Bx()),
+        # Only allocate PML auxilliary fields if needed
+        fCBx = needs_C(sim, Bx()) ? initialize_field_component_array(sim, Bx()) : nothing,
+        fUBx = needs_U(sim, Bx()) ? initialize_field_component_array(sim, Bx()) : nothing,
+        fWBx = needs_W(sim, Bx()) ? initialize_field_component_array(sim, Bx()) : nothing,
+        fHy = initialize_field_component_array(sim, Hy()),
+        fBy = initialize_field_component_array(sim, By()),
+        # Only allocate PML auxilliary fields if needed
+        fCBy = needs_C(sim, By()) ? initialize_field_component_array(sim, By()) : nothing,
+        fUBy = needs_U(sim, By()) ? initialize_field_component_array(sim, By()) : nothing,
+        fWBy = needs_W(sim, By()) ? initialize_field_component_array(sim, By()) : nothing,
+
+        # Only allocate source arrays if we specify sources for that component
+        fSBx = Hx() ∈ sim.source_components ? initialize_field_component_array(sim, Hx()) :
+               nothing,
+        fSBy = Hy() ∈ sim.source_components ? initialize_field_component_array(sim, Hy()) :
+               nothing,
+        fSDz = Ez() ∈ sim.source_components ? initialize_field_component_array(sim, Ez()) :
+               nothing,
     )
 end
 
 """
-    get_fields_args_TM(sim)
+    allocate_fields_TM(sim::SimulationData)
 
-Pull the proper fields arguments for a 2D TM simulation.
+Allocate the proper field components for a 2D TM simulation (where the E-field
+is "in-plane"). This includes Bz/Hz, Dx/Ex, Dy/Ey, and their corresponding
+current sources. 
 """
-function get_fields_args_TM(sim::SimulationData)
+function allocate_fields_TM(sim::SimulationData)
     return (
-        fHz = zeros(get_component_voxel_count(sim, Hz())[1:sim.ndims]...),
-        fBz = zeros(get_component_voxel_count(sim, Bz())[1:sim.ndims]...),
-        fCBz = needs_C(sim, Bz()) ?
-               zeros(get_component_voxel_count(sim, Bz())[1:sim.ndims]...) : nothing,
-        fUBz = needs_U(sim, Bz()) ?
-               zeros(get_component_voxel_count(sim, Bz())[1:sim.ndims]...) : nothing,
-        fWBz = needs_W(sim, Bz()) ?
-               zeros(get_component_voxel_count(sim, Bz())[1:sim.ndims]...) : nothing,
-        fEx = zeros(get_component_voxel_count(sim, Ex())[1:sim.ndims]...),
-        fDx = zeros(get_component_voxel_count(sim, Dx())[1:sim.ndims]...),
-        fCDx = needs_C(sim, Dx()) ?
-               zeros(get_component_voxel_count(sim, Dx())[1:sim.ndims]...) : nothing,
-        fUDx = needs_U(sim, Dx()) ?
-               zeros(get_component_voxel_count(sim, Dx())[1:sim.ndims]...) : nothing,
-        fWDx = needs_W(sim, Dx()) ?
-               zeros(get_component_voxel_count(sim, Dx())[1:sim.ndims]...) : nothing,
-        fEy = zeros(get_component_voxel_count(sim, Ey())[1:sim.ndims]...),
-        fDy = zeros(get_component_voxel_count(sim, Dy())[1:sim.ndims]...),
-        fCDy = needs_C(sim, Dy()) ?
-               zeros(get_component_voxel_count(sim, Dy())[1:sim.ndims]...) : nothing,
-        fUDy = needs_U(sim, Dy()) ?
-               zeros(get_component_voxel_count(sim, Dy())[1:sim.ndims]...) : nothing,
-        fWDy = needs_W(sim, Dy()) ?
-               zeros(get_component_voxel_count(sim, Dy())[1:sim.ndims]...) : nothing,
-        fSDx = Ex() ∈ sim.source_components ?
-               zeros(get_component_voxel_count(sim, Ex())[1:sim.ndims]...) : nothing,
-        fSDy = Ey() ∈ sim.source_components ?
-               zeros(get_component_voxel_count(sim, Ey())[1:sim.ndims]...) : nothing,
-        fSBz = Hz() ∈ sim.source_components ?
-               zeros(get_component_voxel_count(sim, Hz())[1:sim.ndims]...) : nothing,
+        fHz = initialize_field_component_array(sim, Hz()),
+        fBz = initialize_field_component_array(sim, Bz()),
+        # Only allocate PML auxilliary fields if needed
+        fCBz = needs_C(sim, Bz()) ? initialize_field_component_array(sim, Bz()) : nothing,
+        fUBz = needs_U(sim, Bz()) ? initialize_field_component_array(sim, Bz()) : nothing,
+        fWBz = needs_W(sim, Bz()) ? initialize_field_component_array(sim, Bz()) : nothing,
+        fEx = initialize_field_component_array(sim, Ex()),
+        fDx = initialize_field_component_array(sim, Dx()),
+        # Only allocate PML auxilliary fields if needed
+        fCDx = needs_C(sim, Dx()) ? initialize_field_component_array(sim, Dx()) : nothing,
+        fUDx = needs_U(sim, Dx()) ? initialize_field_component_array(sim, Dx()) : nothing,
+        fWDx = needs_W(sim, Dx()) ? initialize_field_component_array(sim, Dx()) : nothing,
+        fEy = initialize_field_component_array(sim, Ey()),
+        fDy = initialize_field_component_array(sim, Dy()),
+        # Only allocate PML auxilliary fields if needed
+        fCDy = needs_C(sim, Dy()) ? initialize_field_component_array(sim, Dy()) : nothing,
+        fUDy = needs_U(sim, Dy()) ? initialize_field_component_array(sim, Dy()) : nothing,
+        fWDy = needs_W(sim, Dy()) ? initialize_field_component_array(sim, Dy()) : nothing,
+
+        # Only allocate source arrays if we specify sources for that component
+        fSDx = Ex() ∈ sim.source_components ? initialize_field_component_array(sim, Ex()) :
+               nothing,
+        fSDy = Ey() ∈ sim.source_components ? initialize_field_component_array(sim, Ey()) :
+               nothing,
+        fSBz = Hz() ∈ sim.source_components ? initialize_field_component_array(sim, Hz()) :
+               nothing,
     )
 end
 
-function init_fields(sim::SimulationData, ::Union{Type{TwoD},Type{ThreeD}})
+"""
+    init_fields(sim::SimulationData, ::Union{Type{TwoD},Type{ThreeD}})
+
+Allocate all the necessary field components for a full 3D simulation.
+"""
+function init_fields(sim::SimulationData, ::Union{Type{ThreeD}})
     sim.fields =
-        Fields{backend_array}(; get_fields_args_TE(sim)..., get_fields_args_TM(sim)...)
+        Fields{AbstractArray}(; allocate_fields_TE(sim)..., allocate_fields_TM(sim)...)
 end
 
+# -------------------------------------------------------------------------- #
+# PML auxilliary field utility functions
+# -------------------------------------------------------------------------- #
+
 """
     needs_C(sim::SimulationData, f::Field)
 
@@ -230,6 +282,9 @@ function needs_W(sim::SimulationData, f::Field)
     return !isnothing(get_pml_conductivity_from_field(sim, f, dir))
 end
 
+# -------------------------------------------------------------------------- #
+# Specific usecase field utility functions
+# -------------------------------------------------------------------------- #
 
 """
     get_fields_for_planewave(
@@ -275,28 +330,3 @@ function get_fields_for_planewave(
     end
     return fields
 end
-
-struct FieldProfile{T}
-    Ex::Vector{T}
-    Ey::Vector{T}
-    Ez::Vector{T}
-    Hx::Vector{T}
-    Hy::Vector{T}
-    Hz::Vector{T}
-    gv_Ex::GridVolume
-    gv_Ey::GridVolume
-    gv_Ez::GridVolume
-    gv_Hx::GridVolume
-    gv_Hy::GridVolume
-    gv_Hz::GridVolume
-end
-
-function overlap_integral(
-    field_profile_1::FieldProfile{T},
-    field_profile_2::FieldProfile{T},
-) where {T}
-    # interpolate fields to center of yee cell
-
-    # perform the integral
-
-end

src/Khronos.jl

@@ -8,6 +8,7 @@ using KernelAbstractions
 using Revise
 using Logging
 using LinearAlgebra
+using OffsetArrays
 
 macro status(exs)
     @logmsg(0, exs)