diff --git a/Project.toml b/Project.toml index 1b64ef7..c84023d 100644 --- a/Project.toml +++ b/Project.toml @@ -7,22 +7,23 @@ version = "0.1.0" ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Einsum = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0" GeometryPrimitives = "17051e67-205e-509e-8301-03b320b998e6" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" 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" +VectorModesolver = "8c544392-5fd9-4640-bea4-e12b3d59d9d1" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" [compat] -Makie = "0.19" +CairoMakie = "0.12" julia = "1.8" [extras] diff --git a/examples/dipole.jl b/examples/dipole.jl index 75f49d7..a748214 100644 --- a/examples/dipole.jl +++ b/examples/dipole.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. # # Simulate a dipole in vacuum. diff --git a/examples/periodic_slab.jl b/examples/periodic_slab.jl index 6444aef..ccbe248 100644 --- a/examples/periodic_slab.jl +++ b/examples/periodic_slab.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. # # Using DFT monitors, compute the transmission of a planewave in 3D through a # simple slab structure. Compare the computed response to the analytic response. diff --git a/examples/sphere.jl b/examples/sphere.jl index 9d7b9be..b8dc0bd 100644 --- a/examples/sphere.jl +++ b/examples/sphere.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. # # Simulate the scattering of a planewave off of a conductive sphere. diff --git a/examples/waveguide.jl b/examples/waveguide.jl new file mode 100644 index 0000000..506bc8c --- /dev/null +++ b/examples/waveguide.jl @@ -0,0 +1,72 @@ +# (c) Meta Platforms, Inc. and affiliates. +# +# Excitation of a dielectric waveguide using a mode source. + +import Khronos +using CairoMakie +using GeometryPrimitives + +Khronos.choose_backend(Khronos.CUDADevice(), Float64) + +function build_waveguide_simulation(λ::Number) + + geometry = [ + Khronos.Object( + Cuboid([0.0, 0.0, 0.0], [100.0, 0.5, 0.22]), + Khronos.Material(ε = 3.4^2), + ), # Waveguide - Si + Khronos.Object( + Cuboid([0.0, 0.0, 0.0], [100.0, 100.0, 100.0]), + Khronos.Material(ε = 1.44^2), + ), # Background material -- SiO2 + ] + + sources = [ + Khronos.ModeSource( + time_profile = Khronos.ContinuousWaveSource(fcen = 1.0 / λ), + frequency = 1.0 / λ, + mode_solver_resolution = 50, + mode_index = 1, + center = [0.0, 0.0, 0.0], + size = [0.0, 2.0, 2.0], + solver_tolerance = 1e-6, + geometry = geometry, + ), + ] + + # Build the simulation object, such that it spans a cube 10μm×10μm×10μm. Place PML 1μm thick. + sim = Khronos.Simulation( + cell_size = [4.0, 4.0, 6.0], + cell_center = [0.0, 0.0, 0.0], + resolution = 25, + geometry = geometry, + sources = sources, + boundaries = [[1.0, 1.0], [1.0, 1.0], [1.0, 1.0]], + ) + return sim +end + +λ = 1.55 +sim = build_waveguide_simulation(λ) + +# scene = Khronos.plot2D( +# sim, +# nothing, +# Khronos.Volume([0.0, 0.0, 0.0], [6.0, 2.0, 2.0]); +# plot_geometry=true, +# ) + + +# scene = Khronos.plot_source(sim, sim.sources[1]) +# save("waveguide_source.png", scene) + +t_end = 40.0; +Khronos.run(sim, until = t_end) + +scene = Khronos.plot2D( + sim, + Khronos.Ey(), + Khronos.Volume([0.0, 0.0, 0.0], [6.0, 4.0, 4.0]); + plot_geometry = true, +) +save("waveguide.png", scene) diff --git a/src/Boundaries.jl b/src/Boundaries.jl index 2d4cc91..9643ac9 100644 --- a/src/Boundaries.jl +++ b/src/Boundaries.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. abstract type Boundary end diff --git a/src/DataStructures.jl b/src/DataStructures.jl index bec9f6a..76ad6e6 100644 --- a/src/DataStructures.jl +++ b/src/DataStructures.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. """ diff --git a/src/Fields.jl b/src/Fields.jl index 539b9de..b1701c5 100644 --- a/src/Fields.jl +++ b/src/Fields.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. # -------------------------------------------------------------------------- # # General field array utility functions diff --git a/src/Geometry.jl b/src/Geometry.jl index f584721..0f468bb 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -1,3 +1,9 @@ +# (c) Meta Platforms, Inc. and affiliates. +# +# Contains all relevant functions for manipulating geometry. +# +# NOTE [smartalecH] Don't use Inf to specify bounds in gometric shapes, as it's +# not compatible with `findfirst`. """ Base.findfirst() @@ -285,3 +291,64 @@ get_mat_conductivity_from_field(sim::SimulationData, ::Electric, ::Y) = sim.geometry_data.σDy get_mat_conductivity_from_field(sim::SimulationData, ::Electric, ::Z) = sim.geometry_data.σDz + +# ---------------------------------------------------------- # +# Material functions +# ---------------------------------------------------------- # + +function get_ε_at_frequency(material, frequency) + ε = zeros(ComplexF64, 3, 3) + + # Account for all the DC terms first + ε .+= isnothing(material.ε) ? 0 : material.ε * I(3) + ε[1, 1] += isnothing(material.εx) ? 0 : material.εx + ε[2, 2] += isnothing(material.εy) ? 0 : material.εy + ε[3, 3] += isnothing(material.εz) ? 0 : material.εz + ε[1, 2] += isnothing(material.εxy) ? 0 : material.εxy + ε[2, 1] += isnothing(material.εxy) ? 0 : material.εxy + ε[1, 3] += isnothing(material.εxz) ? 0 : material.εxz + ε[3, 1] += isnothing(material.εxz) ? 0 : material.εxz + ε[2, 3] += isnothing(material.εyz) ? 0 : material.εyz + ε[3, 2] += isnothing(material.εyz) ? 0 : material.εyz + + # Now handle all the conductivities + ω = 2 * π * frequency + if !isnothing(material.σD) + ε = (1 .+ im * material.σD / ω .* I(3)) .* ε + end + if !isnothing(material.σDx) + ε[1, 1] = (1 + im * material.σDx / ω) * ε[1, 1] + end + if !isnothing(material.σDy) + ε[2, 2] = (1 + im * material.σDy / ω) * ε[2, 2] + end + if !isnothing(material.σDx) + ε[3, 3] = (1 + im * material.σDz / ω) * ε[3, 3] + end + + return ε +end + +function fit_complex_material(ε::Number, frequency::Number) + return Material{Float64}(ε = real(ε), σD = 2 * π * frequency * imag(ε) / real(ε)) +end + +""" + transform_material(material::AbstractArray)::AbstractArray + +Transforms a material tensor `material` by a specified `transformation_matrix`. + +More generally, the susceptibilities χ are transformed to MχMᵀ/|det M|, which +corresponds to [transformation +optics](http://math.mit.edu/~stevenj/18.369/coordinate-transform.pdf) for an +arbitrary curvilinear coordinate transformation with Jacobian matrix M. The +absolute value of the determinant is to prevent inadvertent construction of +left-handed materials +""" +function transform_material( + material::AbstractArray, + transformation_matrix::AbstractArray, +)::AbstractArray + return transformation_matrix * material * transformation_matrix' / + det(transformation_matrix) +end diff --git a/src/Khronos.jl b/src/Khronos.jl index 4c622e0..83e7a37 100644 --- a/src/Khronos.jl +++ b/src/Khronos.jl @@ -1,7 +1,8 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. module Khronos +using Einsum using GeometryPrimitives using Parameters using KernelAbstractions @@ -9,6 +10,7 @@ using Revise using Logging using LinearAlgebra using OffsetArrays +import VectorModesolver macro status(exs) @logmsg(0, exs) @@ -23,10 +25,11 @@ include("load_deps.jl") include("DataStructures.jl") include("utils.jl") +include("Geometry.jl") +include("Mode.jl") include("Boundaries.jl") include("Sources/Sources.jl") include("Fields.jl") -include("Geometry.jl") include("DFT.jl") include("Monitors.jl") include("Timestep.jl") diff --git a/src/Mode.jl b/src/Mode.jl new file mode 100644 index 0000000..ffd1c3c --- /dev/null +++ b/src/Mode.jl @@ -0,0 +1,199 @@ +# (c) Meta Platforms, Inc. and affiliates. +# +# Functions involving mode calculations (e.g. modes sources, mode overlaps) + + +""" + get_mode_profiles(; + frequency::Number, + mode_solver_resolution::Int, + mode_index::Int, + center::Vector{<:Real}, + size::Vector{<:Real}, + solver_tolerance::Number, + geometry::Vector{Object}, + boundary_conditions::Tuple{Int,Int,Int,Int} = (1,1,1,1), + )::VectorModesolver.Mode + +Computes mode profiles for a given `geometry` configuration. + +The mode solver handles all solves in the XY plane. In order to take an +arbitrary cross section in 3D (e.g. an XZ or YZ plane), we need to properly +transform the coordinates, materials, _and_ the fields. This is subtle, because +the H-field is a pseudovector, meaning the coordinate handedness is opposite +that of a normal vector system (it uses the left-hand rule). _Luckily_, the +rotation matrices we use are _proper_ rotations (they have a determinant of +1) +so we can use the same transformation for _all_ quantities. + +To do so, we establish the following conventions: +* Forward propagation is in the +X, +Y, or +Z direction (whichever is aligned + with the monitor's normal direction) +* XY monitors are unchanged (easy) +* XZ monitors are to be rotated 90 degrees around the x-axis, such that the + "new" forward direction is +Z. In this case, the orignal x-axis maps to the + new x-axis, and the original z-axis maps to the new y-axis. +* YZ monitors are first rotated 90 degrees around the z-axis, and then 90 + degrees around the x-axis (a repeat of XZ). In this case, the original + y-axis maps to the new x-axis, and the original z-axis maps to the new + y-axis. +""" +function get_mode_profiles(; + frequency::Number, + mode_solver_resolution::Int, + mode_index::Int, + center::Vector{<:Real}, + size::Vector{<:Real}, + solver_tolerance::Number, + geometry::Vector{Object}, + boundary_conditions::Tuple{Int,Int,Int,Int} = (1, 1, 1, 1), +)::VectorModesolver.Mode + # Note that the volume specified must be a plane, so we need 2 nonzero dimensions + @assert count(x -> x != 0, size) >= 2 "The mode region must be a plane (2D)" + + # Extract the proper grid and transformation function depending on the + # configuration of the monitor. + if size[1] == 0 + # YZ plane + idx_x = 2 # Khronos y maps to VectorModesolver x + idx_y = 3 # Khronos z maps to VectorModesolver y + new_dims = (3, 1, 2, 4) + rotate_material = rotate_YZ_to_XZ + rotate_fields = rotate_XY_to_YZ + elseif size[2] == 0 + # XZ plane + idx_x = 1 # Khronos x maps to VectorModesolver x + idx_y = 3 # Khronos z maps to VectorModesolver y + new_dims = (1, 3, 2, 4) + rotate_material = rotate_XZ_to_XY + rotate_fields = rotate_XY_to_XZ + elseif size[3] == 0 + # XY plane + idx_x = 1 # Khronos x maps to VectorModesolver x + idx_y = 2 # Khronos y maps to VectorModesolver y + new_dims = (1, 2, 3, 4) + rotate_material = x -> x + rotate_fields = x -> x + end + + # Build the domain for the mode solver + step_size = 1.0 / mode_solver_resolution + function _custom_range(idx_dim) + return (-size[idx_dim]/2.0+center[idx_dim]):step_size:(size[idx_dim]/2.0+center[idx_dim]) + end + x_range = _custom_range(idx_x) + y_range = _custom_range(idx_y) + + current_point = copy(center) + + # Build a function that queries the geometry + function ε_profile(x::Float64, y::Float64) + # Since the mode solver is some 2D plane, we need to figure out where we + # are oriented and map back to the full 3D space. + current_point[idx_x] = x + current_point[idx_y] = y + + material_idx = Base.findfirst(current_point, geometry) + if isnothing(material_idx) + # no material specified at current point; return vacuum + return (1.0, 0.0, 0.0, 1.0, 1.0) + else + material_tensor = get_ε_at_frequency(geometry[material_idx].material, frequency) + + # Rotate the material tensor appropriately + material_tensor = rotate_material(material_tensor) + + # Pull the supported tensor elements (εxx, εxy, εyx, εyy, εzz) + return ( + material_tensor[1, 1], + material_tensor[1, 2], + material_tensor[2, 1], + material_tensor[2, 2], + material_tensor[3, 3], + ) + end + end + + mode_solver = VectorModesolver.VectorialModesolver( + λ = 1.0 / frequency, + x = x_range, + y = y_range, + ε = ε_profile, + boundary = boundary_conditions, + ) + mode_solve_results = + VectorModesolver.solve(mode_solver, mode_index, solver_tolerance)[mode_index] + + # rotate the fields back appropriately, and swap the coordinate axes + E_fields = + cat(mode_solve_results.Ex, mode_solve_results.Ey, mode_solve_results.Ez, dims = 4) + H_fields = + cat(mode_solve_results.Hx, mode_solve_results.Hy, mode_solve_results.Hz, dims = 4) + E_fields = permutedims(rotate_fields(E_fields), new_dims) + H_fields = permutedims(rotate_fields(H_fields), new_dims) + + # Build a new mode object with the rotated fields + return VectorModesolver.Mode( + λ = mode_solve_results.λ, + neff = mode_solve_results.neff, + x = x_range, + y = y_range, + Ex = E_fields[:, :, :, 1], + Ey = E_fields[:, :, :, 2], + Ez = E_fields[:, :, :, 3], + Hx = H_fields[:, :, :, 1], + Hy = H_fields[:, :, :, 2], + Hz = H_fields[:, :, :, 3], + ) +end + +# ------------------------------------------------------------------- # +# vector field rotation functions +# ------------------------------------------------------------------- # + +""" + _rotate_vector_field(rotation_matrix::AbstractArray, vector_field::AbstractArray)::AbstractArray + +Perform a batch matrix-vector multiplication on the `vector_field` given a `rotation_matrix`. +""" +function _rotate_vector_field( + rotation_matrix::AbstractArray, + vector_field::AbstractArray, +)::AbstractArray + if ndims(vector_field) == 2 + rotated_vector_field = transform_material(vector_field, rotation_matrix) + else + rotated_vector_field = similar(vector_field) + @einsum rotated_vector_field[k, l, m, i] = + rotation_matrix[i, j] * vector_field[k, l, m, j] + end + + return rotated_vector_field +end + +function rotate_XY_to_XZ(vector_field::AbstractArray) + rotation_matrix = [1 0 0; 0 cos(-pi / 2) -sin(-pi / 2); 0 sin(-pi / 2) cos(-pi / 2)] + return _rotate_vector_field(rotation_matrix, vector_field) +end + +function rotate_XZ_to_XY(vector_field::AbstractArray) + rotation_matrix = [1 0 0; 0 cos(-pi / 2) -sin(-pi / 2); 0 sin(-pi / 2) cos(-pi / 2)]' + return _rotate_vector_field(rotation_matrix, vector_field) +end + +function rotate_XZ_to_YZ(vector_field::AbstractArray) + rotation_matrix = [cos(-pi / 2) -sin(-pi / 2) 0; sin(-pi / 2) cos(-pi / 2) 0; 0 0 1] + return _rotate_vector_field(rotation_matrix, vector_field) +end + +function rotate_YZ_to_XZ(vector_field::AbstractArray) + rotation_matrix = [cos(-pi / 2) -sin(-pi / 2) 0; sin(-pi / 2) cos(-pi / 2) 0; 0 0 1]' + return _rotate_vector_field(rotation_matrix, vector_field) +end + +function rotate_XY_to_YZ(vector_field::AbstractArray) + return vector_field |> rotate_XY_to_XZ |> rotate_XZ_to_YZ +end + +function rotate_YZ_to_XY(vector_field::AbstractArray) + return vector_field |> rotate_YZ_to_XZ |> rotate_XZ_to_XY +end diff --git a/src/Simulation.jl b/src/Simulation.jl index 2bd3b8f..cfb0a86 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. export prepare_simulation!, run, run_until, stop_when_dft_decayed, run_benchmark diff --git a/src/Sources/Sources.jl b/src/Sources/Sources.jl index f93919a..569262d 100644 --- a/src/Sources/Sources.jl +++ b/src/Sources/Sources.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. include("./TimeSources.jl") include("./SpatialSources.jl") diff --git a/src/Sources/SpatialSources.jl b/src/Sources/SpatialSources.jl index af631e0..a0dbcb0 100644 --- a/src/Sources/SpatialSources.jl +++ b/src/Sources/SpatialSources.jl @@ -1,8 +1,14 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. # # Here live the various functions and routines needed to implement a source that # is spatially varying. Importantly, all spatially varying sources take as one # of their parameters the corresponding time profile (temporal dependence). + +export UniformSource, EquivalentSource, PlaneWaveSource, GaussianBeamSource, ModeSource + +# ---------------------------------------------------------- # +# Interface functions +# ---------------------------------------------------------- # # # Each source is defined via the following interface functions. Since the # initialization of the source is *not* done in the main hot loop, and very @@ -10,12 +16,6 @@ # dynamic dispatch. So we'll leverage multiple dispatch as much as possible and # worry about type stability once it's a problem. -export UniformSource, EquivalentSource, PlaneWaveSource, GaussianBeamSource - -# ---------------------------------------------------------- # -# Interface functions -# ---------------------------------------------------------- # - """ get_time_profile(source::Source) @@ -82,15 +82,6 @@ end # ---------------------------------------------------------- # # Uniform source # ---------------------------------------------------------- # - -@with_kw struct UniformSourceData <: Source - time_profile::TimeSource - component::Field - amplitude::Number = 1.0 - center::AbstractVector = [0.0, 0.0, 0.0] - size::AbstractVector = [0.0, 0.0, 0.0] -end - """ UniformSource(; time_profile::TimeSource, @@ -100,7 +91,10 @@ end size::AbstractVector = [0.0, 0.0, 0.0], ) -TBW +Generates a spatially uniform source profile for a given `component`, +`amplitude`, and `time_profile`. + +The size and location are determined by `size` and `center` respectively. """ function UniformSource(; time_profile::TimeSource, @@ -118,6 +112,14 @@ function UniformSource(; ) end +@with_kw struct UniformSourceData <: Source + time_profile::TimeSource + component::Field + amplitude::Number = 1.0 + center::AbstractVector = [0.0, 0.0, 0.0] + size::AbstractVector = [0.0, 0.0, 0.0] +end + function get_source_profile(::UniformSourceData, ::Vector{<:Real}, ::Field) return 1.0 end @@ -129,15 +131,6 @@ end # ---------------------------------------------------------- # # Equivalent source # ---------------------------------------------------------- # - -@with_kw struct EquivalentSourceData <: Source - time_profile::TimeSource - center::AbstractVector = [0.0, 0.0, 0.0] - size::AbstractVector = [0.0, 0.0, 0.0] - amplitude::Number = 1.0 - source_profile::Dict{<:Field,<:Function} -end - """ EquivalentSource(; time_profile::TimeSource, @@ -154,7 +147,7 @@ For the injection to work as expected (i.e. to reproduce the required `E` and `H` fields), the field data must decay by the edges of the source plane, or the source plane must span the entire simulation domain and the fields must match the simulation boundary conditions. The equivalent source currents are fully -defined by the field components tangential to the source plane. For e.g. source +defined by the field components tangential to the source plane. For e.g. sources normal along z, the normal components (Ez and Hz) can be provided but will have no effect on the results, and at least one of the tangential components has to be in the dataset, i.e. at least one of Ex, Ey, Hx, and Hy. @@ -180,60 +173,67 @@ function EquivalentSource(; src_vol = Volume(center = center, size = size) normal_vector = get_normal_vector(src_vol) - transverse_fields = plane_normal_direction(src_vol) - Ex = fields[Ex()] - Ey = fields[Ey()] - Ez = fields[Ez()] - Hx = fields[Hx()] - Hy = fields[Hy()] - Hz = fields[Hz()] + Ex_field = fields[Ex()] + Ey_field = fields[Ey()] + Ez_field = fields[Ez()] + Hx_field = fields[Hx()] + Hy_field = fields[Hy()] + Hz_field = fields[Hz()] - current_sources = dict( - # Electric current J = nHat x H + current_sources = Dict( + # Electric current Ĵ = n̂×Ĥ Ex() => gen_interpolator_from_array( - normal_vector[1] * Hz .- normal_vector[2] * Hy, + normal_vector[2] * Hz_field .- normal_vector[3] * Hy_field, src_vol, ), Ey() => gen_interpolator_from_array( - normal_vector[2] * Hx .- normal_vector[0] * Hz, + normal_vector[3] * Hx_field .- normal_vector[1] * Hz_field, src_vol, ), Ez() => gen_interpolator_from_array( - normal_vector[0] * Hy .- normal_vector[1] * Hx, + normal_vector[1] * Hy_field .- normal_vector[2] * Hx_field, src_vol, ), - # Magnetic current K = - nHat x E + # Magnetic current M̂ = -n̂×Ê Hx() => gen_interpolator_from_array( - normal_vector[2] * Ey .- normal_vector[1] * Ez, + normal_vector[3] * Ey_field .- normal_vector[2] * Ez_field, src_vol, ), Hy() => gen_interpolator_from_array( - normal_vector[0] * Ez .- normal_vector[2] * Ex, + normal_vector[1] * Ez_field .- normal_vector[3] * Ex_field, src_vol, ), Hz() => gen_interpolator_from_array( - normal_vector[1] * Ex .- normal_vector[0] * Ey, + normal_vector[2] * Ex_field .- normal_vector[1] * Ey_field, src_vol, ), ) - # Only create current sources we really need - for component in required_field_components - if (iszero(current_sources[component]) || !(component in transverse_fields)) - delete!(current_sources, component) - end - end + transverse_components = + get_plane_transverse_fields(Volume(center = center, size = size)) + # Filter current_sources based on transverse_components + filtered_current_sources = filter(kv -> kv[1] in transverse_components, current_sources) return EquivalentSourceData( time_profile = time_profile, center = center, size = size, amplitude = amplitude, - source_profile = current_sources, + source_profile = filtered_current_sources, + fields = fields, ) end +@with_kw struct EquivalentSourceData <: Source + time_profile::TimeSource + center::AbstractVector = [0.0, 0.0, 0.0] + size::AbstractVector = [0.0, 0.0, 0.0] + amplitude::Number = 1.0 + source_profile::Dict{<:Field,<:Function} + fields::Dict{<:Field,<:AbstractArray} +end + function get_source_components(source::EquivalentSourceData) return collect(keys(source.source_profile)) end @@ -279,21 +279,6 @@ end # ---------------------------------------------------------- # # Planewave source # ---------------------------------------------------------- # - -@with_kw struct PlaneWaveSourceData <: Source - time_profile::TimeSource - center::AbstractVector{<:Real} = [0.0, 0.0, 0.0] - size::AbstractVector{<:Real} = [0.0, 0.0, 0.0] - amplitude::Number = 1.0 - ε::Number = 1.0 - μ::Number = 1.0 - polarization::AbstractVector{<:Number} - k_vector::AbstractVector{<:Real} - k::Number - normal_direction::AbstractVector{<:Real} -end - - """ PlaneWaveSource(; time_profile::TimeSource, @@ -385,6 +370,19 @@ function PlaneWaveSource(; ) end +@with_kw struct PlaneWaveSourceData <: Source + time_profile::TimeSource + center::AbstractVector{<:Real} = [0.0, 0.0, 0.0] + size::AbstractVector{<:Real} = [0.0, 0.0, 0.0] + amplitude::Number = 1.0 + ε::Number = 1.0 + μ::Number = 1.0 + polarization::AbstractVector{<:Number} + k_vector::AbstractVector{<:Real} + k::Number + normal_direction::AbstractVector{<:Real} +end + """ get_planewave_polarization_scaling( source::PlaneWaveSourceData, @@ -457,10 +455,49 @@ end # Mode source # ---------------------------------------------------------- # -#TODO +function ModeSource(; + frequency::Number, + time_profile::TimeSource, + mode_solver_resolution::Int, + mode_index::Int, + center::Vector{<:Real}, + size::Vector{<:Real}, + solver_tolerance::Number, + amplitude::Number = 1.0, + geometry::Vector{Object}, + boundary_conditions = (1, 1, 1, 1), +)::EquivalentSourceData -@with_kw struct ModeSourceData{N<:Number} <: Source - amplitude::N = 1.0 + # Compute the relevant mode profiles + mode_data = get_mode_profiles(; + frequency = frequency, + mode_solver_resolution = mode_solver_resolution, + mode_index = mode_index, + center = center, + size = size, + solver_tolerance = solver_tolerance, + geometry = geometry, + boundary_conditions = boundary_conditions, + ) + + # Prepare all of the fields at the monitor + fields = Dict( + Ex() => mode_data.Ex, + Ey() => mode_data.Ey, + Ez() => mode_data.Ez, + Hx() => mode_data.Hx, + Hy() => mode_data.Hy, + Hz() => mode_data.Hz, + ) + + # Prepare an equivalent source from the fields + return EquivalentSource( + time_profile = time_profile, + fields = fields, + center = center, + size = size, + amplitude = amplitude, + ) end # ---------------------------------------------------------- # @@ -495,7 +532,7 @@ end polarization::Vector{<:Number}, )::GaussianBeamData -TBW +Generates a Gaussian Beam source. """ function GaussianBeamSource(; time_profile::TimeSource, diff --git a/src/Sources/TimeSources.jl b/src/Sources/TimeSources.jl index 8accecb..ed7787b 100644 --- a/src/Sources/TimeSources.jl +++ b/src/Sources/TimeSources.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. export ContinuousWaveSource, GaussianPulseSource diff --git a/src/Timestep.jl b/src/Timestep.jl index 7f2d22b..dddfd39 100644 --- a/src/Timestep.jl +++ b/src/Timestep.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. # # Here lie the core timestepping kernels for the FDTD algorithm. Thanks to # multiple dispatch, we can simply focus on the fundamental (most complicated) diff --git a/src/Visualization.jl b/src/Visualization.jl index ec80ef2..52293cd 100644 --- a/src/Visualization.jl +++ b/src/Visualization.jl @@ -1,43 +1,14 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. +# +# Functions to visualize the simulation domain, spatial/temporal source +# profiles, and monitor data. using CairoMakie using Statistics -Makie.inline!(true) - -function _get_plane_ranges(gv::GridVolume) - x_range = gv.start_idx[1]:gv.end_idx[1] - y_range = gv.start_idx[2]:gv.end_idx[2] - z_range = gv.start_idx[3]:gv.end_idx[3] - return (x_range, y_range, z_range) -end - -function _pull_fields_from_device(sim::SimulationData, component::Field) - current_fields = get_fields_from_component(sim, component) - array_range = get_component_voxel_count(sim, component) - # Index out the ghost cells and collect to the host - current_fields = Base.Array( - collect(current_fields[1:array_range[1], 1:array_range[2], 1:array_range[3]]), - ) -end - -function _pull_field_slices(sim::SimulationData, component::Field, vol::Volume) - current_fields = _pull_fields_from_device(sim, component) - - vol_x = Volume(vol.center, [0, vol.size[2], vol.size[3]]) - vol_y = Volume(vol.center, [vol.size[1], 0, vol.size[3]]) - vol_z = Volume(vol.center, [vol.size[1], vol.size[2], 0]) - - gv_x = GridVolume(sim, vol_x, component) - gv_y = GridVolume(sim, vol_y, component) - gv_z = GridVolume(sim, vol_z, component) - - YZ_slice = mean(current_fields[_get_plane_ranges(gv_x)...], dims = 1)[1, :, :] - XZ_slice = mean(current_fields[_get_plane_ranges(gv_y)...], dims = 2)[:, 1, :] - XY_slice = mean(current_fields[_get_plane_ranges(gv_z)...], dims = 3)[:, :, 1] +export plot2D, plot_monitor, plot_timesource, plot_source - return (XY_slice, XZ_slice, YZ_slice) -end +Makie.inline!(true) """ plot2D( @@ -48,27 +19,43 @@ end symmetric_field_scaling::Bool = true, ) -Plot each of the three 2D cross sections (XY, XZ, YZ) of a simulation (`sim`) as described by `vol`. +Plot each of the three 2D cross sections (XY, XZ, YZ) of a simulation (`sim`) as +described by `vol`. """ function plot2D( sim::SimulationData, - component::Field, + component::Union{Field,Nothing}, vol::Volume; plot_geometry::Bool = false, symmetric_field_scaling::Bool = true, ) prepare_simulation!(sim) + f = Figure() + vol_size = vol.size - # process requested field components - (XY_slice, XZ_slice, YZ_slice) = _pull_field_slices(sim, component, vol) - vmin = min(minimum(YZ_slice), minimum(XZ_slice), minimum(XY_slice)) - vmax = max(maximum(YZ_slice), maximum(XZ_slice), maximum(XY_slice)) - vmax = max(abs(vmin), vmax) - vmin = -vmax - fields_alpha = 1.0 + # Store all the fields plot configurations in a running dict for convenience. + field_args = Dict() + + # Prepare fields cross sections + if !isnothing(component) + # process requested field components + (XY_slice, XZ_slice, YZ_slice) = _pull_field_slices(sim, component, vol) + vmin = min(minimum(YZ_slice), minimum(XZ_slice), minimum(XY_slice)) + vmax = max(maximum(YZ_slice), maximum(XZ_slice), maximum(XY_slice)) + vmax = max(abs(vmin), vmax) + vmin = -vmax + fields_alpha = 1.0 + + field_args[:transparency] = true + field_args[:colormap] = :bluesreds + if symmetric_field_scaling + field_args[:colorrange] = (vmin, vmax) + end + end + # Prepare geometry cross sections if plot_geometry # for now use εx... although we should generalize this (eps_XY_slice, eps_XZ_slice, eps_YZ_slice) = _pull_field_slices(sim, εx(), vol) @@ -77,67 +64,80 @@ function plot2D( fields_alpha = 0.4 end - # if component isa ε - # colormap = :binary - # else: - # colormap = :bluesreds - # end - # bluesreds - - field_args = Dict() - field_args[:transparency] = true - field_args[:colormap] = :bluesreds - if symmetric_field_scaling - field_args[:colorrange] = (vmin, vmax) - end - - f = Figure(size = (1200, 100)) - + # XY cross section ax_xy = Axis(f[1, 1], title = "XY", xlabel = "X", ylabel = "Y", aspect = DataAspect()) if plot_geometry - heatmap!(ax_xy, eps_XY_slice, colormap = :binary, colorrange = (eps_vmin, eps_vmax)) + xs = range(-vol_size[1] / 2, vol_size[1] / 2, length = size(eps_XY_slice)[1]) + ys = range(-vol_size[2] / 2, vol_size[2] / 2, length = size(eps_XY_slice)[2]) + heatmap!( + ax_xy, + xs, + ys, + eps_XY_slice, + colormap = :binary, + colorrange = (eps_vmin, eps_vmax), + ) + end + if !isnothing(component) + xs = range(-vol_size[1] / 2, vol_size[1] / 2, length = size(XY_slice)[1]) + ys = range(-vol_size[2] / 2, vol_size[2] / 2, length = size(XY_slice)[2]) + hm = heatmap!(ax_xy, xs, ys, XY_slice; field_args...) + #Colorbar(f[:, 2], hm) end - xs = range(-vol_size[1] / 2, vol_size[1] / 2, length = size(XY_slice)[1]) - ys = range(-vol_size[2] / 2, vol_size[2] / 2, length = size(XY_slice)[2]) - hm = heatmap!(ax_xy, xs, ys, XY_slice; field_args...) - Colorbar(f[:, 2], hm) - ax_xz = Axis(f[1, 3], title = "XZ", xlabel = "X", ylabel = "Z", aspect = DataAspect()) + # XZ cross section + ax_xz = Axis(f[1, 2], title = "XZ", xlabel = "X", ylabel = "Z", aspect = DataAspect()) if plot_geometry - heatmap!(ax_xz, eps_XZ_slice, colormap = :binary, colorrange = (eps_vmin, eps_vmax)) + xs = range(-vol_size[1] / 2, vol_size[1] / 2, length = size(eps_XZ_slice)[1]) + zs = range(-vol_size[3] / 2, vol_size[3] / 2, length = size(eps_XZ_slice)[2]) + heatmap!( + ax_xz, + xs, + zs, + eps_XZ_slice, + colormap = :binary, + colorrange = (eps_vmin, eps_vmax), + ) + end + if !isnothing(component) + xs = range(-vol_size[1] / 2, vol_size[1] / 2, length = size(XZ_slice)[1]) + zs = range(-vol_size[3] / 2, vol_size[3] / 2, length = size(XZ_slice)[2]) + hm = heatmap!(ax_xz, xs, zs, XZ_slice; field_args...) + #Colorbar(f[:, 4], hm) end - xs = range(-vol_size[1] / 2, vol_size[1] / 2, length = size(XZ_slice)[1]) - zs = range(-vol_size[3] / 2, vol_size[3] / 2, length = size(XZ_slice)[2]) - hm = heatmap!(ax_xz, xs, zs, XZ_slice; field_args...) - Colorbar(f[:, 4], hm) - ax_yz = Axis(f[1, 5], title = "YZ", xlabel = "Y", ylabel = "Z", aspect = DataAspect()) + # YZ cross section + ax_yz = Axis(f[1, 3], title = "YZ", xlabel = "Y", ylabel = "Z", aspect = DataAspect()) if plot_geometry - heatmap!(ax_yz, eps_YZ_slice, colormap = :binary, colorrange = (eps_vmin, eps_vmax)) + ys = range(-vol_size[2] / 2, vol_size[2] / 2, length = size(eps_YZ_slice)[1]) + zs = range(-vol_size[3] / 2, vol_size[3] / 2, length = size(eps_YZ_slice)[2]) + heatmap!( + ax_yz, + ys, + zs, + eps_YZ_slice, + colormap = :binary, + colorrange = (eps_vmin, eps_vmax), + ) end - ys = range(-vol_size[2] / 2, vol_size[2] / 2, length = size(YZ_slice)[1]) - zs = range(-vol_size[3] / 2, vol_size[3] / 2, length = size(YZ_slice)[2]) - hm = heatmap!(ax_yz, ys, zs, YZ_slice; field_args...) - Colorbar(f[:, end+1], hm) + if !isnothing(component) + ys = range(-vol_size[2] / 2, vol_size[2] / 2, length = size(YZ_slice)[1]) + zs = range(-vol_size[3] / 2, vol_size[3] / 2, length = size(YZ_slice)[2]) + hm = heatmap!(ax_yz, ys, zs, YZ_slice; field_args...) + #Colorbar(f[:, end+1], hm) + end + return f end """ plot_monitor(monitor::DFTMonitor, frequency_idx::Int) -TBW +Plot the field response of a `monitor` at a specified `frequency_idx`. """ function plot_monitor(monitor::DFTMonitor, frequency_idx::Int) freq_data = Base.Array(get_dft_fields(monitor)[:, :, :, frequency_idx]) - if size(freq_data)[1] == 1 - freq_data = freq_data[1, :, :] - elseif size(freq_data)[2] == 1 - freq_data = freq_data[:, 1, :] - elseif size(freq_data)[3] == 1 - freq_data = freq_data[:, :, 1] - else - error("Invalid 3D array (size: $(size(freq_data)))") - end + freq_data = _pull_2D_slice_from_3D(freq_data) freq_data = real(freq_data) @@ -163,14 +163,10 @@ function plot_monitor(monitor::DFTMonitor, frequency_idx::Int) return f end -function plot_timesource(sim::SimulationData, source::Source, frequencies) - plot_timesource(sim, get_time_profile(source), frequencies) -end - """ plot_timesource(sim::SimulationData, time_source::TimeSource, frequencies) -TBW +Plot the time and frequency response of a time-dependent source. """ function plot_timesource(sim::SimulationData, time_source::TimeSource, frequencies) prepare_simulation!(sim) @@ -204,6 +200,119 @@ function plot_timesource(sim::SimulationData, time_source::TimeSource, frequenci return f end +""" + plot_timesource(sim::SimulationData, source::Source, frequencies) + +Plot the time and frequency response of a time-dependent source. +""" +function plot_timesource(sim::SimulationData, source::Source, frequencies) + plot_timesource(sim, get_time_profile(source), frequencies) +end + function plot_source(sim::SimulationData, source::Source) + prepare_simulation!(sim) + source_volume = Volume(center = source.center, size = source.size) + transverse_components = get_plane_transverse_fields(source_volume) + + f = Figure() + num_cols = 2 + + for (i, current_field) in enumerate(transverse_components) + row = div(i - 1, num_cols) + 1 + col = (i - 1) % num_cols + 1 + + ax = Axis(f[row, col*2-1]) # Position axis in every second slot for heatmaps + x_range, y_range = _get_slice_axes(source.fields[current_field], source_volume) + field_profile = _pull_2D_slice_from_3D(source.fields[current_field]) + hm = heatmap!(ax, real(field_profile), colormap = :bluesreds) + cb = Colorbar(f[row, col*2], hm, width = Relative(1 / 8)) + end + + return f + +end + +# ---------------------------------------------------- # +# Utility functions +# ---------------------------------------------------- # + +""" + _pull_2D_slice_from_3D(data::AbstractArray)::AbstractArray + +Given a 3D array with one singeleton dimension, indexes out the corresponding +plane (e.g. for plotting). Similar to a `squeeze`, but safer. +""" +function _pull_2D_slice_from_3D(data::AbstractArray)::AbstractArray + if size(data)[1] == 1 + data = data[1, :, :] + elseif size(data)[2] == 1 + data = data[:, 1, :] + elseif size(data)[3] == 1 + data = data[:, :, 1] + else + error("Invalid 3D array (size: $(size(data)))") + end + +end + +function _get_slice_axes(data::AbstractArray, vol::Volume) + if size(data)[1] == 1 + # YZ plane + idx_x = 2 + idx_y = 3 + elseif size(data)[2] == 1 + idx_x = 1 + idx_y = 3 + elseif size(data)[3] == 1 + idx_x = 1 + idx_y = 2 + else + error("Invalid 3D array (size: $(size(data)))") + end + function _custom_range(idx_dim) + return collect( + range( + -vol.size[idx_dim] / 2.0 + vol.center[idx_dim], + vol.size[idx_dim] / 2.0 + vol.center[idx_dim], + length = size(data)[idx_dim], + ), + ) + end + x_range = _custom_range(idx_x) + y_range = _custom_range(idx_y) + return (x_range, y_range) +end + +function _get_plane_ranges(gv::GridVolume) + x_range = gv.start_idx[1]:gv.end_idx[1] + y_range = gv.start_idx[2]:gv.end_idx[2] + z_range = gv.start_idx[3]:gv.end_idx[3] + return (x_range, y_range, z_range) +end + +function _pull_fields_from_device(sim::SimulationData, component::Field) + current_fields = get_fields_from_component(sim, component) + array_range = get_component_voxel_count(sim, component) + # Index out the ghost cells and collect to the host + current_fields = Base.Array( + collect(current_fields[1:array_range[1], 1:array_range[2], 1:array_range[3]]), + ) +end + +function _pull_field_slices(sim::SimulationData, component::Field, vol::Volume) + current_fields = _pull_fields_from_device(sim, component) + + vol_x = Volume(vol.center, [0, vol.size[2], vol.size[3]]) + vol_y = Volume(vol.center, [vol.size[1], 0, vol.size[3]]) + vol_z = Volume(vol.center, [vol.size[1], vol.size[2], 0]) + gv_x = GridVolume(sim, vol_x, component) + gv_y = GridVolume(sim, vol_y, component) + gv_z = GridVolume(sim, vol_z, component) + + YZ_slice = mean(current_fields[_get_plane_ranges(gv_x)...], dims = 1)[1, :, :] + XZ_slice = mean(current_fields[_get_plane_ranges(gv_y)...], dims = 2)[:, 1, :] + XY_slice = mean(current_fields[_get_plane_ranges(gv_z)...], dims = 3)[:, :, 1] + + return (XY_slice, XZ_slice, YZ_slice) end diff --git a/src/load_deps.jl b/src/load_deps.jl index f7d022b..90e9bee 100644 --- a/src/load_deps.jl +++ b/src/load_deps.jl @@ -1,4 +1,6 @@ -# TODO: Generalize to more backends and allow more user control. +# (c) Meta Platforms, Inc. and affiliates. +# +# Handle all of the backend loading. using CUDA using CUDA.CUDAKernels using Metal diff --git a/src/utils.jl b/src/utils.jl index fd79c51..f454e03 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. """ get_direction_index(direction) @@ -14,17 +14,17 @@ get_direction_index(::Z) = 3 For a given component, get the number of grid voxels throughout the domain. """ -get_component_voxel_count(sim::SimulationData, ::Union{Ex,Dx}) = +get_component_voxel_count(sim::SimulationData, ::Union{Ex,Dx,εx}) = [sim.Nx, sim.Ny + 1, sim.Nz + 1] -get_component_voxel_count(sim::SimulationData, ::Union{Ey,Dy}) = +get_component_voxel_count(sim::SimulationData, ::Union{Ey,Dy,εy}) = [sim.Nx + 1, sim.Ny, sim.Nz + 1] -get_component_voxel_count(sim::SimulationData, ::Union{Ez,Dz}) = +get_component_voxel_count(sim::SimulationData, ::Union{Ez,Dz,εz}) = [sim.Nx + 1, sim.Ny + 1, sim.Nz] -get_component_voxel_count(sim::SimulationData, ::Union{Hx,Bx}) = +get_component_voxel_count(sim::SimulationData, ::Union{Hx,Bx,μx}) = [sim.Nx + 1, sim.Ny, sim.Nz] -get_component_voxel_count(sim::SimulationData, ::Union{Hy,By}) = +get_component_voxel_count(sim::SimulationData, ::Union{Hy,By,μy}) = [sim.Nx, sim.Ny + 1, sim.Nz] -get_component_voxel_count(sim::SimulationData, ::Union{Hz,Bz}) = +get_component_voxel_count(sim::SimulationData, ::Union{Hz,Bz,μz}) = [sim.Nx, sim.Ny, sim.Nz + 1] get_component_voxel_count(sim::SimulationData, ::Center) = [sim.Nx, sim.Ny, sim.Nz] @@ -472,20 +472,20 @@ end """ get_neighbors(x::AbstractArray, point::Float64) -TBW +Returns the two array indices closest to a `point` within a 1D array, `x`. """ function get_neighbors(x::AbstractArray, point::Float64) # get closest point - closest_idx = argmin(abs.(x - point)) + closest_idx = argmin(abs.(x .- point)) # determine which side and get the other point if (x[closest_idx] == point) || (closest_idx) == 1 || (closest_idx == length(x)) left_idx = right_idx = closest_idx elseif x[closest_idx] > point - left_idx = x[closest_idx-1] - right_idx = x[closest_idx] + left_idx = closest_idx - 1 + right_idx = closest_idx else - left_idx = x[closest_idx] - right_idx = x[closest_idx+1] + left_idx = closest_idx + right_idx = closest_idx + 1 end return (left_idx, right_idx) @@ -526,6 +526,10 @@ function bilinear_interpolator( # linearly interpolate top x1 top_val = linear_interpolator(x1, y[:, top_idx], point[1]) + if bottom_idx == top_idx + return top_val + end + # linearly interpolate bottom x1 bottom_val = linear_interpolator(x1, y[:, bottom_idx], point[1]) @@ -554,21 +558,31 @@ function trilinear_interpolator( # bilinearly interpolate top (x1,x2) top_val = bilinear_interpolator(x1, x2, y[:, :, top_idx], point[1:2]) + if bottom_idx == top_idx + return top_val + end + # bilinearly interpolate bottom (x1,x2) bottom_val = bilinear_interpolator(x1, x2, y[:, :, bottom_idx], point[1:2]) # manually linearly interpolate x3 Δz = x3[top_idx] - x3[bottom_idx] + bottom_weight = (point[3] - x3[bottom_idx]) / Δz + top_weight = (x3[top_idx] - point[3]) / Δz + return bottom_weight * bottom_val + top_weight * top_val end """ gen_interpolator_from_array(data::AbstractArray, vol::Volume) -TBW +Given a 3D array (`data`) corresponding to a volume (`vol`) generate a trilinear +interpolator funtion. + +Note that one or more axes may be singleton. """ function gen_interpolator_from_array(data::AbstractArray, vol::Volume) @@ -578,15 +592,23 @@ function gen_interpolator_from_array(data::AbstractArray, vol::Volume) end # pull the grid from the array - ranges = [ - collect( - range( - vol.center[k] - vol.size[k] / 2, - vol.center[k] + vol.size[k] / 2, - length = size(data)[k], - ), - ) for k = 1:3 - ] + ranges = [] + for current_axis = 1:3 + if size(data)[current_axis] == 1 + push!(ranges, [vol.center[current_axis]]) + else + push!( + ranges, + collect( + range( + vol.center[current_axis] - vol.size[current_axis] / 2, + vol.center[current_axis] + vol.size[current_axis] / 2, + length = size(data)[current_axis], + ), + ), + ) + end + end function interpolator(point) return trilinear_interpolator(ranges[1], ranges[2], ranges[3], data, point) diff --git a/test/runtests.jl b/test/runtests.jl index 03d6f70..7cd4220 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. using Test @@ -6,6 +6,8 @@ using Test include("test_interpolation.jl") include("test_grid_volume.jl") + include("test_materials.jl") + include("test_mode.jl") include("test_sources.jl") include("test_timestep.jl") include("test_visualization.jl") diff --git a/test/test_interpolation.jl b/test/test_interpolation.jl index 534fd8f..41fe837 100644 --- a/test/test_interpolation.jl +++ b/test/test_interpolation.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. # test the interpolation and restriction functions inside utils.jl. these are # primarily used to create monitor and source regions, where interpolation and @@ -92,3 +92,16 @@ end @test sum(w) ≈ sx * sy atol = 1e-14 end end + +function test_interpolation(Nx, Ny, Nz, size) + data = rand(Nx, Ny, Nz) + vol = Khronos.Volume(center = [0.0, 0.0, 0.0], size = size) + interp = Khronos.gen_interpolator_from_array(data, vol) + @test all(!isnan.(interp([0.0, 0.0, 0.0]))) +end +@testset "Trilinear interpolation" begin + test_interpolation(24, 12, 1, [2.0, 1.0, 0.0]) + test_interpolation(24, 1, 12, [2.0, 0.0, 1.0]) + test_interpolation(1, 24, 12, [1.0, 2.0, 1.0]) + test_interpolation(24, 12, 12, [2.0, 1.0, 1.0]) +end diff --git a/test/test_materials.jl b/test/test_materials.jl new file mode 100644 index 0000000..20a80fa --- /dev/null +++ b/test/test_materials.jl @@ -0,0 +1,35 @@ +# (c) Meta Platforms, Inc. and affiliates. +# +# Test the functionality of all the material functions. + +import Khronos +using Test + +frequency = 1.0 / 1.55 + +@testset "materials" begin + + ε = 12.4 + simple_isotropic = Khronos.Material(ε = ε) + simple_isotropic_tensor = Khronos.get_ε_at_frequency(simple_isotropic, frequency) + @test isapprox(simple_isotropic_tensor[1, 1], ε) + @test isapprox(simple_isotropic_tensor[2, 2], ε) + @test isapprox(simple_isotropic_tensor[2, 2], ε) + + εx = 1.0 + εy = 2.0 + εz = 3.0 + simple_anisotropic = Khronos.Material(εx = εx, εy = εy, εz = εz) + simple_anisotropic_tensor = Khronos.get_ε_at_frequency(simple_anisotropic, frequency) + @test isapprox(simple_anisotropic_tensor[1, 1], εx) + @test isapprox(simple_anisotropic_tensor[2, 2], εy) + @test isapprox(simple_anisotropic_tensor[3, 3], εz) + + complex_ε = 12.0 + im * 0.1 + complex_mat = Khronos.fit_complex_material(complex_ε, frequency) + simple_complex_tensor = Khronos.get_ε_at_frequency(complex_mat, frequency) + @test isapprox(simple_complex_tensor[1, 1], complex_ε) + @test isapprox(simple_complex_tensor[2, 2], complex_ε) + @test isapprox(simple_complex_tensor[2, 2], complex_ε) + +end diff --git a/test/test_mode.jl b/test/test_mode.jl new file mode 100644 index 0000000..6391fbd --- /dev/null +++ b/test/test_mode.jl @@ -0,0 +1,73 @@ +# (c) Meta Platforms, Inc. and affiliates. +# +# Test the functionality of the mode functions + +import Khronos +using Test + + +function _setup_vector_field(vector_field::AbstractArray) + vector_field[:, :, :, 1] .= 0.0 + vector_field[:, :, :, 2] .= 1.0 + vector_field[:, :, :, 3] .= 2.0 + return vector_field +end + +@testset "mode" begin + # Check that things shuffle like we expect + vector_field = rand(10, 20, 13, 3) + vector_field = _setup_vector_field(vector_field) + XZ_field = similar(vector_field) + XZ_field[:, :, :, 1] .= 0.0 + XZ_field[:, :, :, 2] .= 2.0 + XZ_field[:, :, :, 3] .= -1.0 + @test isapprox(XZ_field, Khronos.rotate_XY_to_XZ(vector_field)) + + vector_field = rand(10, 20, 13, 3) + vector_field = _setup_vector_field(vector_field) + + YZ_field = similar(vector_field) + YZ_field[:, :, :, 1] .= 2.0 + YZ_field[:, :, :, 2] .= 0.0 + YZ_field[:, :, :, 3] .= -1.0 + @test isapprox(YZ_field, Khronos.rotate_XY_to_YZ(vector_field)) + + # All of the mode rotation functions should undo each other. + vector_field = rand(10, 20, 13, 3) + @test isapprox( + vector_field, + vector_field |> Khronos.rotate_XY_to_XZ |> Khronos.rotate_XZ_to_XY, + ) + @test isapprox( + vector_field, + vector_field |> Khronos.rotate_XZ_to_XY |> Khronos.rotate_XY_to_XZ, + ) + + @test isapprox( + vector_field, + vector_field |> Khronos.rotate_XZ_to_YZ |> Khronos.rotate_YZ_to_XZ, + ) + @test isapprox( + vector_field, + vector_field |> Khronos.rotate_YZ_to_XZ |> Khronos.rotate_XZ_to_YZ, + ) + + @test isapprox( + vector_field, + vector_field |> Khronos.rotate_XY_to_YZ |> Khronos.rotate_YZ_to_XY, + ) + @test isapprox( + vector_field, + vector_field |> Khronos.rotate_YZ_to_XY |> Khronos.rotate_XY_to_YZ, + ) + + # Some sanity checks. + @test !isapprox( + vector_field, + vector_field |> Khronos.rotate_XZ_to_XY |> Khronos.rotate_YZ_to_XY, + ) + @test !isapprox( + vector_field, + vector_field |> Khronos.rotate_XZ_to_XY |> Khronos.rotate_XZ_to_XY, + ) +end diff --git a/test/test_sources.jl b/test/test_sources.jl index 99cb9d7..9178f9f 100644 --- a/test/test_sources.jl +++ b/test/test_sources.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. # # Test the functionality of all the sources. diff --git a/test/test_visualization.jl b/test/test_visualization.jl index 29a7033..8adf21a 100644 --- a/test/test_visualization.jl +++ b/test/test_visualization.jl @@ -1,4 +1,4 @@ -# (c) Meta Platforms, Inc. and affiliates. Confidential and proprietary. +# (c) Meta Platforms, Inc. and affiliates. # # Test the functionality of all the sources.