Match theoretical model for electro-optic simulation

Show code cell source Hide code cell source

using Gridap
using Gridap.Geometry
using GridapMakie, CairoMakie

using Femwell.Maxwell.Electrostatic
using Femwell.Thermal

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_kX5XrG".
┌ 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

We start with setting up a square domain. For the boundary conditions, we tag the left and the right side of the model. Furthermore, we create a function which returns 1 indipendent of the tag which is the parameter to descrie the constants of the simplified model.

domain = (-1.0, 1.0, -1.0, 1.0)
partition = (20, 20)
model = simplexify(CartesianDiscreteModel(domain, partition))
labels = get_face_labeling(model)
add_tag!(labels, "left", [1, 3, 7])
add_tag!(labels, "right", [2, 4, 8])
tags = get_face_tag(labels, num_cell_dims(model))
Ω = Triangulation(model)
dΩ = Measure(Ω, 1)
τ = CellField(tags, Ω)
constant_21 = tag -> 21
constant_42 = tag -> 42

Electrostatic

The first step ist to calculate the potential (assuming the electrical conductivity to be k=42). For this we solve the electrostatic equation \(Δϕ = 0\) and define the voltage at two oppositing boundaries to 0V at \(x=0\) and 1V at \(x=1\). The theoretical solution of this function is a linear function.

\[ ϕ(x)=x \]

This would mean the average of the potential over the domain should be

\[ \int ϕ dA / \int 1 dA = 0 \]

p0 = compute_potential(constant_42 ∘ τ, Dict("left" => -1.0, "right" => 1.0))
fig, _, plt = plot(Ω, potential(p0), colormap = :cool)
Colorbar(fig[1, 2], plt)
display(fig)

average_potential = ∑(∫(potential(p0))dΩ) / ∑(∫(1)dΩ)
println("The computed value for the average potential is $average_potential")

The current density can be calculated as

\[ i = k \frac{\mathrm{d}ϕ}{\mathrm{d}ϕ} = 42 \]

and thus the averaged current density over the domain is 42.

average_current_density = ∑(∫(current_density(p0))dΩ) / ∑(∫(1)dΩ) println(“The computed value for the average current density is $average_current_density”)

Using this value, we can caluclate the average power density as

\[ p = k i^2 \]

and thus the averaged power density over the domain is also 42.

average_power_density = ∑(∫(power_density(p0))dΩ) / ∑(∫(1)dΩ)
println("The computed value for the average current density is $average_power_density")

Thermal steady state simulation

Now we calculate the thermal steady state based on the previously calculated locally applied power. For this we chose the thermal conductivity to be \(k_{thermal}=21\) and set the boundaries to 0.

T0 = calculate_temperature(constant_21 ∘ τ, power_density(p0), Dict("boundary" => 0.0))
T_average = ∑(∫(temperature(T0))dΩ) / ∑(∫(1)dΩ)
fig, _, plt = plot(Ω, temperature(T0), colormap = :hot)
Colorbar(fig[1, 2], plt)
display(fig)

Show code cell source Hide code cell source

writevtk(
    Ω,
    "results",
    cellfields = [
        "potential" => potential(p0),
        "current" => current_density(p0),
        "temperature" => temperature(T0),
    ],
)

Thermal transient state simulation

For the simulation of the transient starting with the steady state solution we expect the temperature not to change. Also, we don’t expect it to depend on the thermal thermal diffusitivity.

T_transient = calculate_temperature_transient(
    constant_21 ∘ τ,
    constant_42 ∘ τ,
    power_density(p0),
    Dict("boundary" => 0.0),
    temperature(T0),
    1e-5,
    1e-3,
)
sums = [(t, ∑(∫(u)dΩ) / ∑(∫(1)dΩ) / T_average) for (u, t) in T_transient]
display(lines(sums))

T_transient = calculate_temperature_transient(
    constant_21 ∘ τ,
    constant_42 ∘ τ,
    power_density(p0) * 0,
    Dict("boundary" => 0.0),
    temperature(T0),
    1e-5,
    2e-2,
)
sums = [(t, ∑(∫(u)dΩ) / ∑(∫(1)dΩ) / T_average) for (u, t) in T_transient]
display(lines(sums))

T_transient = calculate_temperature_transient(
    constant_21 ∘ τ,
    constant_42 ∘ τ,
    power_density(p0),
    Dict{String,Float64}(),
    temperature(T0) * 0,
    1e-2,
    1.0,
)
sums = [(t, ∑(∫(u)dΩ) / ∑(∫(1)dΩ)) for (u, t) in T_transient]
display(lines(sums))