3D thermal phase shifter#
Show code cell source
using CairoMakie
using Printf
using Gridap
using GridapGmsh
using Gridap.Geometry
using GridapPETSc
using Femwell.Maxwell.Electrostatic
using Femwell.Thermal
dir = @__DIR__
run(`python $dir/heater_3d_mesh.py`)
Show code cell output
The following 1 direct dependency failed to precompile:
GLMakie [e9467ef8-e4e7-5192-8a1a-b1aee30e663a]
Failed to precompile GLMakie [e9467ef8-e4e7-5192-8a1a-b1aee30e663a] to "/home/runner/.julia/compiled/v1.10/GLMakie/jl_0UT2W4".
┌ Warning: OpenGL/GLFW wasn't loaded correctly or couldn't be initialized.
│ This likely means, you're on a headless server without having OpenGL support setup correctly.
│ Have a look at the troubleshooting section in the readme:
│ https://github.com/MakieOrg/Makie.jl/tree/master/GLMakie#troubleshooting-opengl.
└ @ GLMakie ~/.julia/packages/GLMakie/TH3rf/src/gl_backend.jl:4
ERROR: LoadError: InitError: Exception[GLFW.GLFWError(65550, "X11: The DISPLAY environment variable is missing"), ErrorException("glfwInit failed")]
Stacktrace:
[1] __init__()
@ GLFW ~/.julia/packages/GLFW/wmoTL/src/GLFW.jl:35
[2] run_module_init(mod::Module, i::Int64)
@ Base ./loading.jl:1134
[3] register_restored_modules(sv::Core.SimpleVector, pkg::Base.PkgId, path::String)
@ Base ./loading.jl:1122
[4] _include_from_serialized(pkg::Base.PkgId, path::String, ocachepath::String, depmods::Vector{Any})
@ Base ./loading.jl:1067
[5] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String, build_id::UInt128)
@ Base ./loading.jl:1581
[6] _require(pkg::Base.PkgId, env::String)
@ Base ./loading.jl:1938
[7] __require_prelocked(uuidkey::Base.PkgId, env::String)
@ Base ./loading.jl:1812
[8] #invoke_in_world#3
@ ./essentials.jl:926 [inlined]
[9] invoke_in_world
@ ./essentials.jl:923 [inlined]
[10] _require_prelocked(uuidkey::Base.PkgId, env::String)
@ Base ./loading.jl:1803
[11] macro expansion
@ ./loading.jl:1790 [inlined]
[12] macro expansion
@ ./lock.jl:267 [inlined]
[13] __require(into::Module, mod::Symbol)
@ Base ./loading.jl:1753
[14] #invoke_in_world#3
@ ./essentials.jl:926 [inlined]
[15] invoke_in_world
@ ./essentials.jl:923 [inlined]
[16] require(into::Module, mod::Symbol)
@ Base ./loading.jl:1746
[17] top-level scope
@ ~/.julia/packages/GLMakie/TH3rf/src/gl_backend.jl:2
[18] include(mod::Module, _path::String)
@ Base ./Base.jl:495
[19] include(x::String)
@ GLMakie ~/.julia/packages/GLMakie/TH3rf/src/GLMakie.jl:1
[20] top-level scope
@ ~/.julia/packages/GLMakie/TH3rf/src/GLMakie.jl:87
[21] include
@ ./Base.jl:495 [inlined]
[22] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::Nothing)
@ Base ./loading.jl:2222
[23] top-level scope
@ stdin:3
during initialization of module GLFW
in expression starting at /home/runner/.julia/packages/GLMakie/TH3rf/src/gl_backend.jl:1
in expression starting at /home/runner/.julia/packages/GLMakie/TH3rf/src/GLMakie.jl:1
in expression starting at stdin:
Stacktrace:
[1] pkgerror(msg::String)
@ Pkg.Types /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/Types.jl:70
[2] precompile(ctx::Pkg.Types.Context, pkgs::Vector{Pkg.Types.PackageSpec}; internal_call::Bool, strict::Bool, warn_loaded::Bool, already_instantiated::Bool, timing::Bool, _from_loading::Bool, kwargs::@Kwargs{io::IJulia.IJuliaStdio{Base.PipeEndpoint}})
@ Pkg.API /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/API.jl:1659
[3] precompile(pkgs::Vector{Pkg.Types.PackageSpec}; io::IJulia.IJuliaStdio{Base.PipeEndpoint}, kwargs::@Kwargs{_from_loading::Bool})
@ Pkg.API /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/API.jl:159
[4] precompile
@ /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/API.jl:147 [inlined]
[5] #precompile#114
@ /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/API.jl:146 [inlined]
[6] #invokelatest#2
@ ./essentials.jl:894 [inlined]
[7] invokelatest
@ ./essentials.jl:889 [inlined]
[8] _require(pkg::Base.PkgId, env::String)
@ Base ./loading.jl:1963
[9] __require_prelocked(uuidkey::Base.PkgId, env::String)
@ Base ./loading.jl:1812
[10] #invoke_in_world#3
@ ./essentials.jl:926 [inlined]
[11] invoke_in_world
@ ./essentials.jl:923 [inlined]
[12] _require_prelocked(uuidkey::Base.PkgId, env::String)
@ Base ./loading.jl:1803
[13] macro expansion
@ ./loading.jl:1790 [inlined]
[14] macro expansion
@ ./lock.jl:267 [inlined]
[15] __require(into::Module, mod::Symbol)
@ Base ./loading.jl:1753
[16] #invoke_in_world#3
@ ./essentials.jl:926 [inlined]
[17] invoke_in_world
@ ./essentials.jl:923 [inlined]
[18] require(into::Module, mod::Symbol)
@ Base ./loading.jl:1746
Let’s start with defining the constants for the simulation:
electrical_conductivity = [
"core" => 1e-6,
"box" => 1e-13,
"clad" => 1e-13,
"heater" => 2.3e6,
"via2" => 3.77e7,
"metal2" => 3.77e7,
"via1" => 3.77e7,
"metal3" => 3.77e7,
"metal3#e1" => 3.77e7,
"metal3#e2" => 3.77e7,
]
thermal_conductivities = [
"core" => 90,
"box" => 1.38,
"clad" => 1.38,
"heater" => 400,
"via2" => 400,
"metal2" => 400,
"via1" => 400,
"metal3" => 400,
"metal3#e1" => 400,
"metal3#e2" => 400,
]
thermal_diffisitivities = [
"core" => 90 / 711 / 2330,
"box" => 1.38 / 1000 / 2650,
"clad" => 1.38 / 1000 / 2650,
"heater" => 400 / 376 / 8960,
"via2" => 400 / 376 / 8960,
"metal2" => 400 / 376 / 8960,
"via1" => 400 / 376 / 8960,
"metal3" => 400 / 376 / 8960,
"metal3#e1" => 400 / 376 / 8960,
"metal3#e2" => 400 / 376 / 8960,
]
First we load the mesh, get the labels, define a CellField to evaluate the constants on and create functions which work on the tags
model = GmshDiscreteModel("mesh.msh")
Ω = Triangulation(model)
dΩ = Measure(Ω, 1)
labels = get_face_labeling(model)
tags = get_face_tag(labels, num_cell_dims(model))
τ = CellField(tags, Ω)
electrical_conductivity =
Dict(get_tag_from_name(labels, u) => v for (u, v) in electrical_conductivity)
ϵ_electrical_conductivity(tag) = electrical_conductivity[tag]
thermal_conductivities =
Dict(get_tag_from_name(labels, u) => v for (u, v) in thermal_conductivities)
ϵ_conductivities(tag) = thermal_conductivities[tag]
thermal_diffisitivities =
Dict(get_tag_from_name(labels, u) => v for (u, v) in thermal_diffisitivities)
ϵ_diffisitivities(tag) = thermal_diffisitivities[tag]
The next step is to define the boundary conditions, this can be done simply via julia-dicts:
boundary_potentials = Dict(["metal3#e1___clad" => 0.4, "metal3#e2___clad" => 0.0])
boundary_temperatures =
Dict("metal3#e1___None" => 0.0, "metal3#e2___None" => 0.0, "box___None" => 0.0)
Now we’re ready to do the simulations! First we simulate the electrical potential, then we go on with the temperature
p0 = compute_potential(ϵ_electrical_conductivity ∘ τ, boundary_potentials)
voltage = abs(sum(values(boundary_potentials) .* [-1, 1]))
power = abs(sum(∫(power_density(p0))dΩ))
current = power / voltage
println("Voltage: ", @sprintf("%.2f V", voltage))
println("Current: ", @sprintf("%.2f mA", current * 1e3))
println("Power: ", @sprintf("%.2f mW", power * 1e3))
T0 = calculate_temperature(ϵ_conductivities ∘ τ, power_density(p0), boundary_temperatures)
Yay, now we have both fields simulated! Let’s get the average temperature of the silicon waveguide.
Ω_w = Triangulation(model, tags = "core")
dΩ_w = Measure(Ω_w, 1)
silicon_volume = ∑(∫(1)dΩ_w)
println(
"Average Temperature of silicon waveguide: ",
∑(∫(temperature(T0))dΩ_w) / silicon_volume,
" K",
)
And we write the fields to a vtu-file for visualisation using paraview:
Base.zero(::Type{String}) = ""
writevtk(
Ω,
"results",
cellfields = [
"potential" => potential(p0),
"current" => current_density(p0),
"temperature" => temperature(T0),
"heat flux" => heat_flux(T0),
"tags" => τ,
"tag_name" => (tag -> get_tag_name(labels, tag)) ∘ τ,
],
)
Now let’s get to the transient simulation. We’ll simulate the heatup by starting with zero temperature offset and the power density based on the electrical simulation. For the cooldown simulation we start with the steady state temperature profile and no heating.
figure = Figure()
ax = Axis(
figure[1, 1],
ylabel = "Average silicon waveguide temperature / K",
xlabel = "Time / ms",
)
for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", 0, 1)]
uₕₜ = calculate_temperature_transient(
ϵ_conductivities ∘ τ,
ϵ_diffisitivities ∘ τ,
power_density(p0) * power_factor,
boundary_temperatures,
temperature(T0) * temperature_factor,
2e-7,
2e-4,
)
#createpvd("poisson_transient_solution_$label") do pvd
# for (uₕ, t) in uₕₜ
# pvd[t] = createvtk(
# Ω,
# "poisson_transient_solution_$t" * ".vtu",
# cellfields = ["u" => uₕ],
# )
# end
#end
sums = [(t, ∑(∫(u)dΩ_w) / silicon_volume) for (u, t) in uₕₜ]
t, s = getindex.(sums, 1), getindex.(sums, 2)
lines!(ax, t * 1e3, s, label = label)
end
axislegend(position = :rc)
display(figure)
For more complex examples we’d need a iterative solver, i.e.
GridapPETSc.with(args = split("-ksp_type cg -pc_type mg -ksp_monitor")) do
# Here the code
end
And then add to each of the calculate-calls the sovler parameter:
solver = PETScLinearSolver()