Merge branch 'feature/parallelization'

Author: HoBeZwe

.codecov.yml

@@ -3,4 +3,4 @@ coverage:
     project:
       default:
         target: 100%
-        threshold: 1%
+        threshold: 2%

.github/workflows/CI.yml

@@ -4,9 +4,12 @@ on:
     branches:
       - 'master'
       - 'release-'
+      - 'feature/*'
     tags: 
       - '*'
   pull_request:
+    branches:
+      - '*'
   release:
     types: [published]
 

Project.toml

@@ -1,12 +1,14 @@
 name = "SphericalScattering"
 uuid = "1a9ea918-b599-4f1f-bd9a-d681e8bb5b3e"
 authors = ["Bernd Hofmann <[email protected]> and contributors"]
-version = "0.7.2"
+version = "0.8.0"
 
 [deps]
 LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
 PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
@@ -19,7 +21,9 @@ SphericalScatteringExt = "PlotlyJS"
 [compat]
 LegendrePolynomials = "0.3, 0.4"
 LinearAlgebra = "1"
+OhMyThreads = "0.8.3"
 PlotlyJS = "0.18"
+ProgressMeter = "1.10.4"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"
 julia = "1.10"

src/SphericalScattering.jl

@@ -20,6 +20,8 @@ const ε0 = 8.8541878176e-12  # default permittivity
 using SpecialFunctions, LegendrePolynomials
 using LinearAlgebra
 using StaticArrays
+using OhMyThreads
+using ProgressMeter
 
 
 

src/UniformField/scattered.jl

@@ -8,10 +8,15 @@ function scatteredfield(sphere::Sphere, excitation::UniformField, quantity::Fiel
 
     F = zeros(fieldType(quantity), size(quantity.locations))
 
+    points = quantity.locations
+    p = progress(length(points))
+
     # --- compute field in Cartesian representation
-    for (ind, point) in enumerate(quantity.locations)
-        F[ind] = scatteredfield(sphere, excitation, point, quantity; parameter=parameter)
+    @tasks for ind in eachindex(points)
+        F[ind] = scatteredfield(sphere, excitation, points[ind], quantity; parameter=parameter)
+        next!(p)
     end
+    finish!(p)
 
     return F
 end

src/dataHandling.jl

@@ -38,3 +38,8 @@ struct Parameter
 end
 
 Parameter() = Parameter(-1, 1e-12)
+
+# global setting for the style of the progress bar
+function progress(numIter::Int)
+    return Progress(numIter; barglyphs=BarGlyphs("[=> ]"), color=:white)
+end

src/dipoles/incident.jl

@@ -35,8 +35,8 @@ function field(excitation::Dipole, point, quantity::ElectricField; parameter::Pa
 
     Il = excitation.amplitude
     k  = wavenumber(excitation)
-    ε  = excitation.embedding.ε
-    μ  = excitation.embedding.μ
+    ε = excitation.embedding.ε
+    μ = excitation.embedding.μ
 
     r0 = excitation.position
     p  = excitation.orientation
@@ -87,8 +87,8 @@ function field(excitation::HertzianDipole, point, quantity::FarField; parameter:
 
     Il = excitation.amplitude
     k  = wavenumber(excitation)
-    ε  = excitation.embedding.ε
-    μ  = excitation.embedding.μ
+    ε = excitation.embedding.ε
+    μ = excitation.embedding.μ
 
     r0 = excitation.position
     p  = excitation.orientation

src/planeWave/scattered.jl

@@ -12,10 +12,14 @@ function scatteredfield(sphere::Sphere, excitation::PlaneWave, quantity::Field;
     # --- rotate coordinates
     points = rotate(excitation, quantity.locations; inverse=true)
 
+    p = progress(length(points))
+
     # --- compute field in Cartesian representation
-    for (ind, point) in enumerate(points)
-        F[ind] = scatteredfield(sphere, excitation, point, quantity; parameter=parameter)
+    @tasks for ind in eachindex(points)
+        F[ind] = scatteredfield(sphere, excitation, points[ind], quantity; parameter=parameter)
+        next!(p)
     end
+    finish!(p)
 
     # --- rotate resulting field
     rotate!(excitation, F; inverse=false)
@@ -42,10 +46,9 @@ function scatteredfield(sphere::Sphere, excitation::PlaneWave, point, quantity::
 
     eps = parameter.relativeAccuracy
 
-    ST = SVector{3,Complex{T}}
-    F = ST(0.0, 0.0, 0.0)
+    F = SVector{3,Complex{T}}(0.0, 0.0, 0.0)
 
-    A₀ = amplitude(sphere, excitation::PlaneWave, quantity, r)
+    A₀ = amplitude(sphere, excitation, quantity, r)
 
     A₀ == 0.0 && return F # Inside of PEC return zero field
 
@@ -80,7 +83,30 @@ function scatteredfield(sphere::Sphere, excitation::PlaneWave, point, quantity::
 
     end
 
-    return convertSpherical2Cartesian(A₀ .* F, point_sph)
+    Fin = inside(sphere, excitation, point, quantity; parameter=parameter)
+
+    return convertSpherical2Cartesian(A₀ .* F, point_sph) + Fin
+end
+
+
+function inside(sphere::Sphere, excitation::PlaneWave{T,R,C}, point, quantity::Field; parameter) where {T,R,C}
+
+    return SVector{3,Complex{R}}(0.0, 0.0, 0.0) # no correction needed
+end
+
+function inside(sphere::DielectricSphere, excitation::PlaneWave{T,R,C}, point, quantity::FarField; parameter) where {T,R,C}
+
+    return SVector{3,Complex{R}}(0.0, 0.0, 0.0) # no correction needed
+end
+
+function inside(sphere::DielectricSphere, excitation::PlaneWave{T,R,C}, point, quantity::Field; parameter) where {T,R,C}
+
+    # inside the sphere the incident field has to be substracted to get only the scattered part
+    if norm(point) < sphere.radius
+        return -field(excitation, point, quantity; parameter=parameter)
+    end
+
+    return SVector{3,Complex{R}}(0.0, 0.0, 0.0)
 end
 
 

src/ringCurrent/incident.jl

@@ -122,10 +122,10 @@ function field(excitation::RingCurrent, point, quantity::MagneticField; paramete
 
     eps = parameter.relativeAccuracy
 
-    Hr  = Complex{T}(0.0) # initialize
+    Hr   = Complex{T}(0.0) # initialize
     Hϑ  = Complex{T}(0.0) # initialize
     δHr = T(Inf)
-    n   = -1
+    n    = -1
 
     r = point_sph[1]
 

src/ringCurrent/scattered.jl

@@ -19,10 +19,14 @@ function scatteredfield(sphere::PECSphere, excitation::RingCurrent, quantity::Fi
     # --- translate/rotate coordinates
     points = rotate(excitation, quantity.locations; inverse=true)
 
+    p = progress(length(points))
+
     # --- compute field in Cartesian representation
-    for (ind, point) in enumerate(points)
-        F[ind] = scatteredfield(sphere, exc, point, fieldType; parameter=parameter)
+    @tasks for ind in eachindex(points)
+        F[ind] = scatteredfield(sphere, exc, points[ind], fieldType; parameter=parameter)
+        next!(p)
     end
+    finish!(p)
 
     # --- rotate resulting field
     rotate!(excitation, F; inverse=false)

src/sphere.jl

@@ -210,7 +210,7 @@ end
 wavenumber(sp::Sphere, ex::Excitation, r)
 
 Returns the wavenumber at radius `r` in the sphere `sp`.
-If this part is PEC, k=0 is rweturned.
+If this part is PEC, k=0 is returned.
 """
 function wavenumber(sp::PECSphere, ex::Excitation, r)
     ε = ex.embedding.ε

src/sphericalModes/scattered.jl

@@ -19,10 +19,14 @@ function scatteredfield(sphere::PECSphere, excitation::SphericalMode, quantity::
 
     γ = scatterCoeff(sphere, excitation, excitation.n, wavenumber(excitation) * sphere.radius)
 
+    p = progress(length(points))
+
     # --- compute field in Cartesian representation
-    for (ind, point) in enumerate(points)
-        F[ind] = γ * scatteredfield(sphere, exc, point, fieldType; parameter=parameter)
+    @tasks for ind in eachindex(points)
+        F[ind] = γ * scatteredfield(sphere, exc, points[ind], fieldType; parameter=parameter)
+        next!(p)
     end
+    finish!(p)
 
     # --- rotate resulting field
     # rotate!(F, excitation.rotation)

src/utils.jl

@@ -4,13 +4,13 @@
 
 Convenience function returning points (in Cartesian and spherical coordinates) on a spherical grid at a distance 'r' and a resolution in degrees. 
 """
-function sphericalGridPoints(; r=1.0, resolution=5)
+function sphericalGridPoints(; r=1.0, resolution=5, center=SVector(0.0, 0.0, 0.0))
 
     ϑ = range(0.0; stop=π, length=round(Int, 180 / resolution))  # default: 5° steps
     ϕ = range(0.0; stop=2π, length=round(Int, 360 / resolution)) # default: 5° steps 
 
-    points_cart = [SVector(r * cos(φ) * sin(θ), r * sin(φ) * sin(θ), r * cos(θ)) for θ in ϑ, φ in ϕ]
-    points_sph  = [SVector(r, θ, φ) for θ in ϑ, φ in ϕ]
+    points_cart = [SVector(r * cos(φ) * sin(θ), r * sin(φ) * sin(θ), r * cos(θ)) + center for θ in ϑ, φ in ϕ]
+    points_sph  = [SVector(r, θ, φ) + center for θ in ϑ, φ in ϕ]
 
     return points_cart, points_sph
 end

test/planeWave_dielectric.jl

@@ -77,7 +77,7 @@
 
     # E-Field
     EF₂ = scatteredfield(sp, ex, ElectricField(points_cartNF))
-    EF₁ = scatteredfield(sp, ex, ElectricField(points_cartNF_inside))
+    EF₁ = field(sp, ex, ElectricField(points_cartNF_inside))
 
     diff_EF₂ = norm.(EF₂ - EF₂MoM) ./ maximum(norm.(EF₂))  # worst case error
     diff_EF₁ = norm.(EF₁ - EF₁MoM) ./ maximum(norm.(EF₁))  # worst case error
@@ -90,7 +90,7 @@
     HF₁MoM = hfield(𝓣k1, -(1 / η1)^2 .* m, RT, 𝓚k1, -j, RT, points_cartNF_inside)
 
     HF₂ = scatteredfield(sp, ex, MagneticField(points_cartNF))
-    HF₁ = scatteredfield(sp, ex, MagneticField(points_cartNF_inside))
+    HF₁ = field(sp, ex, MagneticField(points_cartNF_inside))
 
     diff_HF₂ = norm.(HF₂ - HF₂MoM) ./ maximum(norm.(HF₂))  # worst case error
     diff_HF₁ = norm.(HF₁ - HF₁MoM) ./ maximum(norm.(HF₁))  # worst case error