add julia mode solver to benchmark

docs/benchmarks/mode_solver.py

@@ -27,6 +27,7 @@
 import pandas as pd
 import shapely
 import shapely.affinity
+from juliacall import Main as jl
 from skfem import Basis, ElementTriP0
 from skfem.io.meshio import from_meshio
 
@@ -43,6 +44,7 @@
 neff_values_paper = [3.412022, 3.412126, 3.412279, 3.412492, 3.412774, 3.413132, 3.413571, 3.414100]
 neff_values_femwell_slepc = []
 neff_values_femwell_scipy = []
+neff_values_julia = []
 
 for slab_thickness in slab_thicknesses:
     if slab_thickness == 0.0:
@@ -81,10 +83,37 @@
             polygons, resolutions, default_resolution_max=0.3, filename="mesh.msh"
         )
     )
+    refractive_indices = {"core": 3.44, "box": 3.40, "clad": 1}
+    jl.refractive_indices = refractive_indices
+
+    julia_script = [
+        "using Pkg",
+        'Pkg.activate("../../")',
+        "using Gridap",
+        "using Gridap.Geometry",
+        "using GridapGmsh",
+        "using Femwell.Maxwell.Waveguide",
+        "using GridapMakie, CairoMakie",
+        "import PhysicalConstants.CODATA2018: c_0, μ_0, ε_0",
+        "import Unitful: ustrip",
+        'model = GmshDiscreteModel("mesh.msh")',
+        "Ω = Triangulation(model)",
+        "labels = get_face_labeling(model)",
+        "τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)",
+        'refractive_indices = Dict("core" => 3.44^2, "box" => 3.40^2, "clad" => 1^2)',
+        "ε(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in refractive_indices)[tag]",
+        "modes = calculate_modes(model, ε ∘ τ,  λ = 1.15, order = 1)",
+        "neffs = [real(n_eff(mode)) for mode in modes]",
+    ]
+
+    for line in julia_script:
+        jl.seval(line)
+
+    neff_values_julia.append(jl.neffs[0])
 
     basis0 = Basis(mesh, ElementTriP0())
     epsilon = basis0.zeros()
-    for subdomain, n in {"core": 3.44, "box": 3.40, "clad": 1}.items():
+    for subdomain, n in refractive_indices.items():
         epsilon[basis0.get_dofs(elements=subdomain)] = n**2
 
     modes = compute_modes(basis0, epsilon, wavelength=1.15, num_modes=1, order=2, solver="slepc")
@@ -93,6 +122,7 @@
     modes = compute_modes(basis0, epsilon, wavelength=1.15, num_modes=1, order=2, solver="scipy")
     neff_values_femwell_scipy.append(np.real(modes[0].n_eff))
 
+
 pd.DataFrame(
     {
         "slab_thickness": slab_thicknesses,
@@ -105,13 +135,17 @@
         "difference scipy": (
             f"{n1-n2:.6f}" for n1, n2 in zip(neff_values_paper, neff_values_femwell_scipy)
         ),
+        "calculated value julia": (f"{n:.6f}" for n in neff_values_julia),
+        "difference julia": (
+            f"{n1-n2:.6f}" for n1, n2 in zip(neff_values_paper, neff_values_julia)
+        ),
     }
 ).style.apply(
     lambda differences: [
         "background-color: green" if abs(float(difference)) < 5e-6 else "background-color: red"
         for difference in differences
     ],
-    subset=["difference slepc", "difference scipy"],
+    subset=["difference slepc", "difference scipy", "difference julia"],
 )
 # -