From 084c5b7660c78c83860986e8550ad1a073bcdce2 Mon Sep 17 00:00:00 2001 From: Ian Hammond Date: Sat, 27 Jan 2024 16:05:57 -0500 Subject: [PATCH 1/8] heater_3d mesh and fields plots --- docs/julia/heater_3d.jl | 92 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 90 insertions(+), 2 deletions(-) diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl index b51087ba..534d1fca 100644 --- a/docs/julia/heater_3d.jl +++ b/docs/julia/heater_3d.jl @@ -19,6 +19,7 @@ # %% tags=["remove-stderr", "hide-input", "hide-output"] using CairoMakie +using GridapMakie using Printf using Gridap @@ -32,6 +33,46 @@ using Femwell.Thermal dir = @__DIR__ run(`python $dir/heater_3d_mesh.py`) +# %% [markdown] +# We can visualize the mesh using GridapMakie: + +# %% tags=["remove-stderr", "hide-input", "hide-output"] +function visualize_mesh(model::DiscreteModel) + volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] + colors = ["grey", "black","orange","green", "green", "brown", "brown", "blue", "blue"] + alphas = [1.0, 1.0, 0.2, 0.5, 0.5, 0.1, 0.1, 0.5, 0.5] + fig = Figure(fontsize=20) + xlims = (-20e-6, 60e-6) + ylims = (-5e-6, 5e-6) + zlims = (-5e-6, 5e-6) + aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) + tickfunc = values->(x->string(x*1e6)).(values) + + # Mesh + ∂Ω = BoundaryTriangulation(model, tags=volume_tags) + ax = fig[1,2] = Axis3( + fig[1,2], aspect=(aspect, 1.0, 1.0), + xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)\n\n", title="Heater Mesh", + limits=(xlims, ylims, zlims), titlesize=28, + xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, + azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 + ) + plt = plot!(fig[1,2], ∂Ω, shading=true, title="Mesh", fontsize=20) + for (i, tag) in enumerate(volume_tags) + ∂Ω_region = BoundaryTriangulation(model, tags=[tag]) + wireframe!(∂Ω_region, color=(colors[i], alphas[i]), linewidth=0.1, fontsize=20) + end + wf = [PolyElement(color=color) for color in colors] + Legend(fig[1,1], wf, volume_tags, halign=:left, valign=:top, framevisible=false, labelsize=20) + + return fig +end + +# %% tags=["remove-stderr"] +model = GmshDiscreteModel("mesh.msh") +fig = visualize_mesh(model) +display(fig) + # %% [markdown] # Let's start with defining the constants for the simulation: @@ -77,7 +118,6 @@ thermal_diffisitivities = [ # and create functions which work on the tags # %% tags=["remove-stderr", "hide-output"] -model = GmshDiscreteModel("mesh.msh") Ω = Triangulation(model) dΩ = Measure(Ω, 1) @@ -135,7 +175,55 @@ println( ) # %% [markdown] -# And we write the fields to a file for visualisation using paraview: +# And we plot the fields using GridapMakie: + +# %% tags=["remove-stderr", "hide-input", "hide-output"] +function plot_fields(model, potential, current_density, temperature, heat_flux) + + # Set up figure + volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] + fontsize = 10 + fig = Figure(fontsize=fontsize, resolution=(800, 500)) + xlims = (-20e-6, 60e-6) + ylims = (-5e-6, 5e-6) + zlims = (-5e-6, 5e-6) + aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) + tickfunc = values->(x->string(x*1e6)).(values) + + # Set up four quadrants + g11 = fig[1, 1] = GridLayout() + g12 = fig[1, 2] = GridLayout() + g21 = fig[2, 1] = GridLayout() + g22 = fig[2, 2] = GridLayout() + + # Plot fields + gs = [g11, g12, g21, g22] + fs = [potential, current_density, temperature, heat_flux] + labels = ["Potential", "Current", "Temperature", "Heat Flux"] + cmaps = [:viridis, :viridis, Reverse(:RdBu), Reverse(:RdBu)] + cranges = [(0, 0.5), (0, 1.75e10), (0, 120.0), (3.0e4, 4e9)] + Ω = Triangulation(model, tags=volume_tags) + for (g, f, label, cmap, crange) ∈ zip(gs, fs, labels, cmaps, cranges) + Axis3( + g[1,1], aspect=(aspect, 1.0, 1.0), + xlabel="x (μm)", ylabel="y (μm)", zlabel="z (μm)", title=label, + limits=(xlims, ylims, zlims), titlesize=18, + xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, + azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 + ) + plt = plot!(g[1,1], Ω, f, shading=true, colorrange=crange, colormap=cmap, fontsize=fontsize) + Colorbar(g[1,2], plt, label=label, vertical=true, flipaxis = false) + colsize!(g, 1, Aspect(1, 1.25)) + end + return fig +end + +# %% tags=["remove-stderr"] +fig = plot_fields(model, potential(p0), current_density(p0), temperature(T0), heat_flux(T0)) +display(fig) + +# %% [markdown] +# Alternatively, we write the fields to a file for visualisation using paraview: # %% tags=["remove-stderr", "hide-output"] writevtk( From bacaa97e2e843e530e3ebcdbc63eb6ea2153c7e8 Mon Sep 17 00:00:00 2001 From: Ian Hammond Date: Sat, 27 Jan 2024 19:33:50 -0500 Subject: [PATCH 2/8] transient animation --- docs/julia/heater_3d.jl | 96 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 87 insertions(+), 9 deletions(-) diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl index 534d1fca..909f5b8b 100644 --- a/docs/julia/heater_3d.jl +++ b/docs/julia/heater_3d.jl @@ -69,7 +69,7 @@ function visualize_mesh(model::DiscreteModel) end # %% tags=["remove-stderr"] -model = GmshDiscreteModel("mesh.msh") +model = GmshDiscreteModel("$dir/mesh.msh") fig = visualize_mesh(model) display(fig) @@ -182,8 +182,8 @@ function plot_fields(model, potential, current_density, temperature, heat_flux) # Set up figure volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] - fontsize = 10 - fig = Figure(fontsize=fontsize, resolution=(800, 500)) + fontsize = 10*2 + fig = Figure(fontsize=fontsize, resolution=(800*2, 500*2)) xlims = (-20e-6, 60e-6) ylims = (-5e-6, 5e-6) zlims = (-5e-6, 5e-6) @@ -206,8 +206,8 @@ function plot_fields(model, potential, current_density, temperature, heat_flux) for (g, f, label, cmap, crange) ∈ zip(gs, fs, labels, cmaps, cranges) Axis3( g[1,1], aspect=(aspect, 1.0, 1.0), - xlabel="x (μm)", ylabel="y (μm)", zlabel="z (μm)", title=label, - limits=(xlims, ylims, zlims), titlesize=18, + xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)\n", title=label, + limits=(xlims, ylims, zlims), titlesize=18*2, xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 ) @@ -228,7 +228,7 @@ display(fig) # %% tags=["remove-stderr", "hide-output"] writevtk( Ω, - "results", + "$dir/results", cellfields = [ "potential" => potential(p0), "current" => current_density(p0), @@ -252,6 +252,69 @@ ax = Axis( xlabel = "Time / ms", ) +# %% [markdown] +# We can plot the transient +# %% tags=["remove-stderr", "hide-input", "hide-output"] +function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; label="") + # Set up figure + volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] + fig = Figure() + xlims = (-20e-6, 60e-6) + ylims = (-5e-6, 5e-6) + zlims = (-5e-6, 5e-6) + aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) + tickfunc = values->(x->string(x*1e6)).(values) + Ω = Triangulation(model, tags=volume_tags) + fig = Figure() + + # Set up axis + it = iterate(uₕₜ) + uh = Any[nothing] + time = Any[nothing] + state = Any[nothing] + (uh[1], time[1]), state[1] = it + t_uh_time = Observable(time[1]) + timetitle = lift(time -> label * @sprintf(" t = %.2f ms", time * 1e3), t_uh_time) + Axis3( + fig[1,1], aspect=(aspect, 1.0, 1.0), + xlabel="x (μm)", ylabel="y (μm)", zlabel="z (μm)", title=timetitle, + limits=(xlims, ylims, zlims), + xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, + azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 + ) + + # Plot new iteration + u = lift(t_uh_time) do t + while (time[1] < t) || (time[1] ≈ t) + if time[1] ≈ t + @show time[1] + return uh[1] + end + it = iterate(uₕₜ, state[1]) + if isnothing(it) + return uh[1] + end + (uh[1], time[1]), state[1] = it + end + return uh[1] + end + + # Create the actual plot + plt = plot!(fig[1,1], Ω, u, shading=true, colorrange=(0, 110.0), colormap=Reverse(:RdBu)) + Colorbar(fig[1,2], plt, vertical=true, label="Temperature (K)") + colsize!(fig.layout, 1, Aspect(1, 1.0)) + t_uh_time, fig +end + +# %% [markdown] +# Here we set the parameters for the transient simulation +# %% tags=["remove-stderr", "hide-output"] +Δt = 2e-7 +t_end = 2e-4 +plot_every = 200 +timestamps = collect(Δt:Δt*plot_every:t_end) +total_render_time = 5.0 +fps = ceil(Int, length(timestamps) / total_render_time) for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", 0, 1)] uₕₜ = calculate_temperature_transient( ϵ_conductivities ∘ τ, @@ -259,10 +322,17 @@ for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", power_density(p0) * power_factor, boundary_temperatures, temperature(T0) * temperature_factor, - 2e-7, - 2e-4, + Δt, + t_end, ) + # Create animation + t_uh_time, fig = plot_time(uₕₜ; label) + record(fig, "$dir/animation-$label.gif", timestamps, framerate=fps) do this_t + t_uh_time[] = this_t + end + + # Optionally write the solution to a file #createpvd("poisson_transient_solution_$label") do pvd # for (uₕ, t) in uₕₜ # pvd[t] = createvtk( @@ -278,6 +348,14 @@ for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", lines!(ax, t * 1e3, s, label = label) end +# %% [markdown] +# If we display the animations we can visualize the heatup and cooldown of the device: +# ![heatup animation](animation-heatup.gif) +# ![cooldown animation](animation-cooldown.gif) + +# %% [markdown] +# Finally, we can plot an average of the temperature over time: +# %% tags=["remove-stderr"] axislegend(position = :rc) display(figure) @@ -291,7 +369,7 @@ display(figure) # end # ``` # -# And then add to each of the calculate-calls the sovler parameter: +# And then add to each of the calculate-calls the solver parameter: # # ```julia # solver = PETScLinearSolver() From 9d19d63c9016416e497be2aa7ca8beb03cbff2cc Mon Sep 17 00:00:00 2001 From: Ian Hammond Date: Sun, 28 Jan 2024 12:14:40 -0500 Subject: [PATCH 3/8] removed show time in transient plot function --- docs/julia/heater_3d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl index 909f5b8b..2f07370d 100644 --- a/docs/julia/heater_3d.jl +++ b/docs/julia/heater_3d.jl @@ -287,7 +287,6 @@ function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; la u = lift(t_uh_time) do t while (time[1] < t) || (time[1] ≈ t) if time[1] ≈ t - @show time[1] return uh[1] end it = iterate(uₕₜ, state[1]) From 0a67789b4655a748974e5a36f44da270acd5b4b2 Mon Sep 17 00:00:00 2001 From: Ian Hammond Date: Sun, 28 Jan 2024 12:29:29 -0500 Subject: [PATCH 4/8] added some spacing for some plot labels --- docs/julia/heater_3d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl index 2f07370d..61399ffc 100644 --- a/docs/julia/heater_3d.jl +++ b/docs/julia/heater_3d.jl @@ -206,7 +206,7 @@ function plot_fields(model, potential, current_density, temperature, heat_flux) for (g, f, label, cmap, crange) ∈ zip(gs, fs, labels, cmaps, cranges) Axis3( g[1,1], aspect=(aspect, 1.0, 1.0), - xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)\n", title=label, + xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)\n\n", title=label, limits=(xlims, ylims, zlims), titlesize=18*2, xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 @@ -277,7 +277,7 @@ function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; la timetitle = lift(time -> label * @sprintf(" t = %.2f ms", time * 1e3), t_uh_time) Axis3( fig[1,1], aspect=(aspect, 1.0, 1.0), - xlabel="x (μm)", ylabel="y (μm)", zlabel="z (μm)", title=timetitle, + xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)", title=timetitle, limits=(xlims, ylims, zlims), xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 From 5bf31df962b99c62d7a4414b58380ab119ae61d5 Mon Sep 17 00:00:00 2001 From: Ian Hammond Date: Sat, 17 Feb 2024 17:16:24 -0500 Subject: [PATCH 5/8] interactive plots --- Project.toml | 5 + docs/julia/heater_3d.jl | 260 ++++++++++++++++++++++------------------ 2 files changed, 151 insertions(+), 114 deletions(-) diff --git a/Project.toml b/Project.toml index 2016bdc2..1cd7edfb 100644 --- a/Project.toml +++ b/Project.toml @@ -5,21 +5,26 @@ version = "0.0.0-DEV" [deps] Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" +Bonito = "824d6782-a2ef-11e9-3a09-e5662e0c26f8" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8" GridapMakie = "41f30b06-6382-4b60-a5f7-79d86b35bf5d" GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1" +Gumbo = "708ec375-b3d6-5a57-a7ce-8257bf98657a" +HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3" IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MUMPS = "55d2b088-9f4e-11e9-26c0-150b02ea6a46" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" Pickle = "fbb45041-c46e-462f-888f-7c521cafbc2c" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" [compat] Arpack = "=0.5.3" diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl index 61399ffc..828ed4a9 100644 --- a/docs/julia/heater_3d.jl +++ b/docs/julia/heater_3d.jl @@ -18,7 +18,6 @@ # # 3D thermal phase shifter # %% tags=["remove-stderr", "hide-input", "hide-output"] -using CairoMakie using GridapMakie using Printf @@ -30,18 +29,57 @@ using GridapPETSc using Femwell.Maxwell.Electrostatic using Femwell.Thermal -dir = @__DIR__ -run(`python $dir/heater_3d_mesh.py`) +using IJulia +using CairoMakie +using WGLMakie +using Gumbo, HTTP +using Bonito, Markdown +WGLMakie.activate!() -# %% [markdown] -# We can visualize the mesh using GridapMakie: +Makie.inline!(true) # Make sure to inline plots into Documenter output! + +dir = @__DIR__ +# run(`python $dir/heater_3d_mesh.py`) # %% tags=["remove-stderr", "hide-input", "hide-output"] -function visualize_mesh(model::DiscreteModel) +js_counter = 1 +function separate_javascript!(element::HTMLElement) + for (i, child) in enumerate(element.children) + if typeof(child) <: HTMLElement{:script} + src_raw = replace(child.attributes["src"], "data:application/javascript;base64," => "") + js_content = HTTP.base64decode(src_raw) + global js_counter += 1 + js_filename = "./js/js-$(js_counter).js" + js_write = "$dir/../_build/html/julia/js/js-$(js_counter).js" + write(js_write, js_content) + element.children[i] = HTMLElement{:script}(Vector{HTMLNode}[], element, Dict{String,String}("src"=>js_filename, "type"=>"module")) + end + if typeof(child) <: HTMLElement + separate_javascript!(child) + end + end + return element +end + +function separate_javascript(html_source::AbstractString) + if !isdir("$dir/../_build/html/") + mkdir("$dir/../_build/html/") + end + if !isdir("$dir/../_build/html/julia/") + mkdir("$dir/../_build/html/julia/") + end + if !isdir("$dir/../_build/html/julia/js/") + mkdir("$dir/../_build/html/julia/js/") + end + parsed_html = parsehtml(html_source) + div_element = separate_javascript!(parsed_html.root.children[2].children[1]) +end + +function visualize_mesh(model) volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] - colors = ["grey", "black","orange","green", "green", "brown", "brown", "blue", "blue"] - alphas = [1.0, 1.0, 0.2, 0.5, 0.5, 0.1, 0.1, 0.5, 0.5] - fig = Figure(fontsize=20) + colors = [:grey, :black, :orange, :green, :green, :brown, :brown, :blue, :blue] + alphas = [1.0, 1.0, 0.2, 0.5, 0.5, 0.1, 0.1, 0.75, 0.75] + fig = Figure(fontsize=16) xlims = (-20e-6, 60e-6) ylims = (-5e-6, 5e-6) zlims = (-5e-6, 5e-6) @@ -49,7 +87,6 @@ function visualize_mesh(model::DiscreteModel) tickfunc = values->(x->string(x*1e6)).(values) # Mesh - ∂Ω = BoundaryTriangulation(model, tags=volume_tags) ax = fig[1,2] = Axis3( fig[1,2], aspect=(aspect, 1.0, 1.0), xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)\n\n", title="Heater Mesh", @@ -57,21 +94,102 @@ function visualize_mesh(model::DiscreteModel) xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 ) - plt = plot!(fig[1,2], ∂Ω, shading=true, title="Mesh", fontsize=20) + camera = Camera3D(ax.scene) + for (i, tag) in enumerate(volume_tags) ∂Ω_region = BoundaryTriangulation(model, tags=[tag]) - wireframe!(∂Ω_region, color=(colors[i], alphas[i]), linewidth=0.1, fontsize=20) + wireframe!(fig[1,2], ∂Ω_region, color=colors[i], alpha=alphas[i], linewidth=0.1, fontsize=16) end wf = [PolyElement(color=color) for color in colors] - Legend(fig[1,1], wf, volume_tags, halign=:left, valign=:top, framevisible=false, labelsize=20) + Legend(fig[1,1], wf, volume_tags, halign=:left, valign=:top, framevisible=false, labelsize=16) + + html_str = string(DOM.div(fig)) + html_source_new = separate_javascript(html_str) + display("text/html", string(html_source_new)) +end - return fig +function visualize_field(model, field, label; plotargs...) + volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] + fontsize = 10 + fig = Figure(fontsize=fontsize, resolution=(800, 500)) + xlims = (-20e-6, 60e-6) + ylims = (-5e-6, 5e-6) + zlims = (-5e-6, 5e-6) + aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) + tickfunc = values->(x->string(x*1e6)).(values) + + # Setup plot + Ω = Triangulation(model, tags=volume_tags) + ax = fig[1,1] = Axis3( + fig[1,1], aspect=(aspect, 1.0, 1.0), + xlabel="x (μm)", ylabel="y (μm)", zlabel="z (μm)", title=label, + limits=(xlims, ylims, zlims), titlesize=18, + xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, + azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 + ) + camera = Camera3D(ax.scene) + + # Plot + plt = plot!(fig[1,1], Ω, field, shading=true; plotargs...) + Colorbar(fig[1,2], plt, label=label, vertical=true, flipaxis = false) + + html_str = string(DOM.div(fig)) + html_source_new = separate_javascript(html_str) + display("text/html", string(html_source_new)) +end + +function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; label="") + # Set up figure + volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] + fig = CairoMakie.Figure() + xlims = (-20e-6, 60e-6) + ylims = (-5e-6, 5e-6) + zlims = (-5e-6, 5e-6) + aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) + tickfunc = values->(x->string(x*1e6)).(values) + Ω = Triangulation(model, tags=volume_tags) + + # Set up axis + it = iterate(uₕₜ) + uh = Any[nothing] + time = Any[nothing] + state = Any[nothing] + (uh[1], time[1]), state[1] = it + t_uh_time = Observable(time[1]) + timetitle = lift(time -> label * @sprintf(" t = %.2f ms", time * 1e3), t_uh_time) + Axis3( + fig[1,1], aspect=(aspect, 1.0, 1.0), + xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)", title=timetitle, + limits=(xlims, ylims, zlims), + xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, + azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 + ) + + # Plot new iteration + u = lift(t_uh_time) do t + while (time[1] < t) || (time[1] ≈ t) + if time[1] ≈ t + return uh[1] + end + it = iterate(uₕₜ, state[1]) + if isnothing(it) + return uh[1] + end + (uh[1], time[1]), state[1] = it + end + return uh[1] + end + + # Create the actual plot + plt = plot!(fig[1,1], Ω, u, colorrange=(0, 110.0), colormap=Reverse(:RdBu)) + Colorbar(fig[1,2], plt, vertical=true, label="Temperature (K)") + colsize!(fig.layout, 1, Aspect(1, 1.0)) + t_uh_time, fig end # %% tags=["remove-stderr"] model = GmshDiscreteModel("$dir/mesh.msh") -fig = visualize_mesh(model) -display(fig) +visualize_mesh(model) # %% [markdown] # Let's start with defining the constants for the simulation: @@ -174,53 +292,17 @@ println( " K", ) -# %% [markdown] -# And we plot the fields using GridapMakie: - -# %% tags=["remove-stderr", "hide-input", "hide-output"] -function plot_fields(model, potential, current_density, temperature, heat_flux) +# %% tags=["remove-stderr"] +visualize_field(model, potential(p0), "Potential"; colormap=:viridis, colorrange=(0, 0.5), fontsize=14) - # Set up figure - volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] - fontsize = 10*2 - fig = Figure(fontsize=fontsize, resolution=(800*2, 500*2)) - xlims = (-20e-6, 60e-6) - ylims = (-5e-6, 5e-6) - zlims = (-5e-6, 5e-6) - aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) - tickfunc = values->(x->string(x*1e6)).(values) +# %% tags=["remove-stderr"] +visualize_field(model, current_density(p0), "Current Density"; colormap=:viridis, colorrange=(0, 1.75e10), fontsize=14) - # Set up four quadrants - g11 = fig[1, 1] = GridLayout() - g12 = fig[1, 2] = GridLayout() - g21 = fig[2, 1] = GridLayout() - g22 = fig[2, 2] = GridLayout() - - # Plot fields - gs = [g11, g12, g21, g22] - fs = [potential, current_density, temperature, heat_flux] - labels = ["Potential", "Current", "Temperature", "Heat Flux"] - cmaps = [:viridis, :viridis, Reverse(:RdBu), Reverse(:RdBu)] - cranges = [(0, 0.5), (0, 1.75e10), (0, 120.0), (3.0e4, 4e9)] - Ω = Triangulation(model, tags=volume_tags) - for (g, f, label, cmap, crange) ∈ zip(gs, fs, labels, cmaps, cranges) - Axis3( - g[1,1], aspect=(aspect, 1.0, 1.0), - xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)\n\n", title=label, - limits=(xlims, ylims, zlims), titlesize=18*2, - xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, - azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 - ) - plt = plot!(g[1,1], Ω, f, shading=true, colorrange=crange, colormap=cmap, fontsize=fontsize) - Colorbar(g[1,2], plt, label=label, vertical=true, flipaxis = false) - colsize!(g, 1, Aspect(1, 1.25)) - end - return fig -end +# %% tags=["remove-stderr"] +visualize_field(model, temperature(T0), "Temperature"; colormap=Reverse(:RdBu), colorrange=(0, 120.0), fontsize=14) # %% tags=["remove-stderr"] -fig = plot_fields(model, potential(p0), current_density(p0), temperature(T0), heat_flux(T0)) -display(fig) +visualize_field(model, heat_flux(T0), "Heat Flux"; colormap=Reverse(:RdBu), colorrange=(3.0e4, 4e9), fontsize=14) # %% [markdown] # Alternatively, we write the fields to a file for visualisation using paraview: @@ -245,6 +327,7 @@ writevtk( # For the cooldown simulation we start with the steady state temperature profile and no heating. # %% tags=["remove-stderr"] +CairoMakie.activate!() figure = Figure() ax = Axis( figure[1, 1], @@ -252,59 +335,6 @@ ax = Axis( xlabel = "Time / ms", ) -# %% [markdown] -# We can plot the transient -# %% tags=["remove-stderr", "hide-input", "hide-output"] -function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; label="") - # Set up figure - volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] - fig = Figure() - xlims = (-20e-6, 60e-6) - ylims = (-5e-6, 5e-6) - zlims = (-5e-6, 5e-6) - aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) - tickfunc = values->(x->string(x*1e6)).(values) - Ω = Triangulation(model, tags=volume_tags) - fig = Figure() - - # Set up axis - it = iterate(uₕₜ) - uh = Any[nothing] - time = Any[nothing] - state = Any[nothing] - (uh[1], time[1]), state[1] = it - t_uh_time = Observable(time[1]) - timetitle = lift(time -> label * @sprintf(" t = %.2f ms", time * 1e3), t_uh_time) - Axis3( - fig[1,1], aspect=(aspect, 1.0, 1.0), - xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)", title=timetitle, - limits=(xlims, ylims, zlims), - xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, - azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 - ) - - # Plot new iteration - u = lift(t_uh_time) do t - while (time[1] < t) || (time[1] ≈ t) - if time[1] ≈ t - return uh[1] - end - it = iterate(uₕₜ, state[1]) - if isnothing(it) - return uh[1] - end - (uh[1], time[1]), state[1] = it - end - return uh[1] - end - - # Create the actual plot - plt = plot!(fig[1,1], Ω, u, shading=true, colorrange=(0, 110.0), colormap=Reverse(:RdBu)) - Colorbar(fig[1,2], plt, vertical=true, label="Temperature (K)") - colsize!(fig.layout, 1, Aspect(1, 1.0)) - t_uh_time, fig -end - # %% [markdown] # Here we set the parameters for the transient simulation # %% tags=["remove-stderr", "hide-output"] @@ -314,6 +344,8 @@ plot_every = 200 timestamps = collect(Δt:Δt*plot_every:t_end) total_render_time = 5.0 fps = ceil(Int, length(timestamps) / total_render_time) +html_sources = String[] +lins = CairoMakie.Lines[] for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", 0, 1)] uₕₜ = calculate_temperature_transient( ϵ_conductivities ∘ τ, @@ -344,7 +376,7 @@ for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", 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) + push!(lins,lines!(ax, t * 1e3, s, label = label)) end # %% [markdown] @@ -355,7 +387,7 @@ end # %% [markdown] # Finally, we can plot an average of the temperature over time: # %% tags=["remove-stderr"] -axislegend(position = :rc) +Legend(figure[1, 1], lins, ["heatup","cooldown"], halign=:right, valign=:center, labelsize=16) display(figure) From 572107e4305e5aa11fa5bfff4ac397202ffc5fa3 Mon Sep 17 00:00:00 2001 From: Ian Hammond Date: Thu, 22 Feb 2024 13:23:22 -0500 Subject: [PATCH 6/8] separate js code & display methods --- docs/julia/heater_3d.jl | 254 ++++++++++++++++++++++++++-------------- docs/julia/js.jl | 77 ++++++++++++ 2 files changed, 240 insertions(+), 91 deletions(-) create mode 100644 docs/julia/js.jl diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl index 828ed4a9..dcdb8d0b 100644 --- a/docs/julia/heater_3d.jl +++ b/docs/julia/heater_3d.jl @@ -21,6 +21,9 @@ using GridapMakie using Printf +include("./js.jl") +using .JS: FigureContainer + using Gridap using GridapGmsh using Gridap.Geometry @@ -31,9 +34,9 @@ using Femwell.Thermal using IJulia using CairoMakie +using Markdown + using WGLMakie -using Gumbo, HTTP -using Bonito, Markdown WGLMakie.activate!() Makie.inline!(true) # Make sure to inline plots into Documenter output! @@ -42,112 +45,139 @@ dir = @__DIR__ # run(`python $dir/heater_3d_mesh.py`) # %% tags=["remove-stderr", "hide-input", "hide-output"] -js_counter = 1 -function separate_javascript!(element::HTMLElement) - for (i, child) in enumerate(element.children) - if typeof(child) <: HTMLElement{:script} - src_raw = replace(child.attributes["src"], "data:application/javascript;base64," => "") - js_content = HTTP.base64decode(src_raw) - global js_counter += 1 - js_filename = "./js/js-$(js_counter).js" - js_write = "$dir/../_build/html/julia/js/js-$(js_counter).js" - write(js_write, js_content) - element.children[i] = HTMLElement{:script}(Vector{HTMLNode}[], element, Dict{String,String}("src"=>js_filename, "type"=>"module")) - end - if typeof(child) <: HTMLElement - separate_javascript!(child) - end - end - return element -end - -function separate_javascript(html_source::AbstractString) - if !isdir("$dir/../_build/html/") - mkdir("$dir/../_build/html/") - end - if !isdir("$dir/../_build/html/julia/") - mkdir("$dir/../_build/html/julia/") - end - if !isdir("$dir/../_build/html/julia/js/") - mkdir("$dir/../_build/html/julia/js/") - end - parsed_html = parsehtml(html_source) - div_element = separate_javascript!(parsed_html.root.children[2].children[1]) -end function visualize_mesh(model) - volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] + volume_tags = [ + "box", + "core", + "heater", + "via2", + "via1", + "metal3", + "metal2", + "metal3#e1", + "metal3#e2", + ] colors = [:grey, :black, :orange, :green, :green, :brown, :brown, :blue, :blue] alphas = [1.0, 1.0, 0.2, 0.5, 0.5, 0.1, 0.1, 0.75, 0.75] - fig = Figure(fontsize=16) + fig = Figure(fontsize = 16) xlims = (-20e-6, 60e-6) ylims = (-5e-6, 5e-6) zlims = (-5e-6, 5e-6) - aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) - tickfunc = values->(x->string(x*1e6)).(values) + aspect = (xlims[2] - xlims[1]) / (ylims[2] - ylims[1]) + tickfunc = values -> (x -> string(x * 1e6)).(values) # Mesh - ax = fig[1,2] = Axis3( - fig[1,2], aspect=(aspect, 1.0, 1.0), - xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)\n\n", title="Heater Mesh", - limits=(xlims, ylims, zlims), titlesize=28, - xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, - azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 - ) - camera = Camera3D(ax.scene) + ax = + fig[1, 2] = Axis3( + fig[1, 2], + aspect = (aspect, 1.0, 1.0), + xlabel = "x (μm)", + ylabel = "\ny (μm)", + zlabel = "z (μm)\n\n", + title = "Heater Mesh", + limits = (xlims, ylims, zlims), + titlesize = 28, + xtickformat = tickfunc, + ytickformat = tickfunc, + ztickformat = tickfunc, + azimuth = 51 / 40 * pi + pi / 20, + elevation = pi / 8 + pi / 10, + ) + Camera3D(ax.scene) for (i, tag) in enumerate(volume_tags) - ∂Ω_region = BoundaryTriangulation(model, tags=[tag]) - wireframe!(fig[1,2], ∂Ω_region, color=colors[i], alpha=alphas[i], linewidth=0.1, fontsize=16) + ∂Ω_region = BoundaryTriangulation(model, tags = [tag]) + wireframe!( + fig[1, 2], + ∂Ω_region, + color = colors[i], + alpha = alphas[i], + linewidth = 0.1, + fontsize = 16, + ) end - wf = [PolyElement(color=color) for color in colors] - Legend(fig[1,1], wf, volume_tags, halign=:left, valign=:top, framevisible=false, labelsize=16) + wf = [PolyElement(color = color) for color in colors] + Legend( + fig[1, 1], + wf, + volume_tags, + halign = :left, + valign = :top, + framevisible = false, + labelsize = 16, + ) - html_str = string(DOM.div(fig)) - html_source_new = separate_javascript(html_str) - display("text/html", string(html_source_new)) + display(FigureContainer(fig)) end function visualize_field(model, field, label; plotargs...) - volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] + volume_tags = [ + "box", + "core", + "heater", + "via2", + "via1", + "metal3", + "metal2", + "metal3#e1", + "metal3#e2", + ] fontsize = 10 - fig = Figure(fontsize=fontsize, resolution=(800, 500)) + fig = Figure(fontsize = fontsize, resolution = (800, 500)) xlims = (-20e-6, 60e-6) ylims = (-5e-6, 5e-6) zlims = (-5e-6, 5e-6) - aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) - tickfunc = values->(x->string(x*1e6)).(values) + aspect = (xlims[2] - xlims[1]) / (ylims[2] - ylims[1]) + tickfunc = values -> (x -> string(x * 1e6)).(values) # Setup plot - Ω = Triangulation(model, tags=volume_tags) - ax = fig[1,1] = Axis3( - fig[1,1], aspect=(aspect, 1.0, 1.0), - xlabel="x (μm)", ylabel="y (μm)", zlabel="z (μm)", title=label, - limits=(xlims, ylims, zlims), titlesize=18, - xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, - azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 - ) - camera = Camera3D(ax.scene) + Ω = Triangulation(model, tags = volume_tags) + ax = + fig[1, 1] = Axis3( + fig[1, 1], + aspect = (aspect, 1.0, 1.0), + xlabel = "x (μm)", + ylabel = "y (μm)", + zlabel = "z (μm)", + title = label, + limits = (xlims, ylims, zlims), + titlesize = 18, + xtickformat = tickfunc, + ytickformat = tickfunc, + ztickformat = tickfunc, + azimuth = 51 / 40 * pi + pi / 20, + elevation = pi / 8 + pi / 10, + ) + Camera3D(ax.scene) # Plot - plt = plot!(fig[1,1], Ω, field, shading=true; plotargs...) - Colorbar(fig[1,2], plt, label=label, vertical=true, flipaxis = false) - - html_str = string(DOM.div(fig)) - html_source_new = separate_javascript(html_str) - display("text/html", string(html_source_new)) + plt = plot!(fig[1, 1], Ω, field, shading = true; plotargs...) + Colorbar(fig[1, 2], plt, label = label, vertical = true, flipaxis = false) + + display(FigureContainer(fig)) end -function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; label="") +function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; label = "") # Set up figure - volume_tags = ["box","core","heater","via2","via1","metal3","metal2","metal3#e1","metal3#e2"] + volume_tags = [ + "box", + "core", + "heater", + "via2", + "via1", + "metal3", + "metal2", + "metal3#e1", + "metal3#e2", + ] fig = CairoMakie.Figure() xlims = (-20e-6, 60e-6) ylims = (-5e-6, 5e-6) zlims = (-5e-6, 5e-6) - aspect = (xlims[2]-xlims[1])/(ylims[2]-ylims[1]) - tickfunc = values->(x->string(x*1e6)).(values) - Ω = Triangulation(model, tags=volume_tags) + aspect = (xlims[2] - xlims[1]) / (ylims[2] - ylims[1]) + tickfunc = values -> (x -> string(x * 1e6)).(values) + Ω = Triangulation(model, tags = volume_tags) # Set up axis it = iterate(uₕₜ) @@ -158,11 +188,18 @@ function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; la t_uh_time = Observable(time[1]) timetitle = lift(time -> label * @sprintf(" t = %.2f ms", time * 1e3), t_uh_time) Axis3( - fig[1,1], aspect=(aspect, 1.0, 1.0), - xlabel="x (μm)", ylabel="\ny (μm)", zlabel="z (μm)", title=timetitle, - limits=(xlims, ylims, zlims), - xtickformat=tickfunc, ytickformat=tickfunc, ztickformat=tickfunc, - azimuth=51/40 * pi + pi/20, elevation=pi/8 + pi / 10 + fig[1, 1], + aspect = (aspect, 1.0, 1.0), + xlabel = "x (μm)", + ylabel = "\ny (μm)", + zlabel = "z (μm)", + title = timetitle, + limits = (xlims, ylims, zlims), + xtickformat = tickfunc, + ytickformat = tickfunc, + ztickformat = tickfunc, + azimuth = 51 / 40 * pi + pi / 20, + elevation = pi / 8 + pi / 10, ) # Plot new iteration @@ -181,8 +218,8 @@ function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; la end # Create the actual plot - plt = plot!(fig[1,1], Ω, u, colorrange=(0, 110.0), colormap=Reverse(:RdBu)) - Colorbar(fig[1,2], plt, vertical=true, label="Temperature (K)") + plt = plot!(fig[1, 1], Ω, u, colorrange = (0, 110.0), colormap = Reverse(:RdBu)) + Colorbar(fig[1, 2], plt, vertical = true, label = "Temperature (K)") colsize!(fig.layout, 1, Aspect(1, 1.0)) t_uh_time, fig end @@ -293,16 +330,44 @@ println( ) # %% tags=["remove-stderr"] -visualize_field(model, potential(p0), "Potential"; colormap=:viridis, colorrange=(0, 0.5), fontsize=14) +visualize_field( + model, + potential(p0), + "Potential"; + colormap = :viridis, + colorrange = (0, 0.5), + fontsize = 14, +) # %% tags=["remove-stderr"] -visualize_field(model, current_density(p0), "Current Density"; colormap=:viridis, colorrange=(0, 1.75e10), fontsize=14) +visualize_field( + model, + current_density(p0), + "Current Density"; + colormap = :viridis, + colorrange = (0, 1.75e10), + fontsize = 14, +) # %% tags=["remove-stderr"] -visualize_field(model, temperature(T0), "Temperature"; colormap=Reverse(:RdBu), colorrange=(0, 120.0), fontsize=14) +visualize_field( + model, + temperature(T0), + "Temperature"; + colormap = Reverse(:RdBu), + colorrange = (0, 120.0), + fontsize = 14, +) # %% tags=["remove-stderr"] -visualize_field(model, heat_flux(T0), "Heat Flux"; colormap=Reverse(:RdBu), colorrange=(3.0e4, 4e9), fontsize=14) +visualize_field( + model, + heat_flux(T0), + "Heat Flux"; + colormap = Reverse(:RdBu), + colorrange = (3.0e4, 4e9), + fontsize = 14, +) # %% [markdown] # Alternatively, we write the fields to a file for visualisation using paraview: @@ -359,7 +424,7 @@ for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", # Create animation t_uh_time, fig = plot_time(uₕₜ; label) - record(fig, "$dir/animation-$label.gif", timestamps, framerate=fps) do this_t + record(fig, "$dir/animation-$label.gif", timestamps, framerate = fps) do this_t t_uh_time[] = this_t end @@ -376,7 +441,7 @@ for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", sums = [(t, ∑(∫(u)dΩ_w) / silicon_volume) for (u, t) in uₕₜ] t, s = getindex.(sums, 1), getindex.(sums, 2) - push!(lins,lines!(ax, t * 1e3, s, label = label)) + push!(lins, lines!(ax, t * 1e3, s, label = label)) end # %% [markdown] @@ -387,7 +452,14 @@ end # %% [markdown] # Finally, we can plot an average of the temperature over time: # %% tags=["remove-stderr"] -Legend(figure[1, 1], lins, ["heatup","cooldown"], halign=:right, valign=:center, labelsize=16) +Legend( + figure[1, 1], + lins, + ["heatup", "cooldown"], + halign = :right, + valign = :center, + labelsize = 16, +) display(figure) diff --git a/docs/julia/js.jl b/docs/julia/js.jl new file mode 100644 index 00000000..7b0bf675 --- /dev/null +++ b/docs/julia/js.jl @@ -0,0 +1,77 @@ +module JS + +using Gumbo, HTTP, Bonito +using Makie, WGLMakie +export FigureContainer + +abstract type BackendType end +struct WGLMakieType <: BackendType end + +function backend_type(::Module) + if Makie.current_backend() === WGLMakie + return WGLMakieType + else + return BackendType + end +end + +struct FigureContainer{T<:BackendType} + fig::Figure +end +FigureContainer(fig::Figure) = FigureContainer{backend_type(Makie.current_backend())}(fig) + +js_counter = 1 +dir = @__DIR__ + +function separate_javascript!(element::HTMLElement) + for (i, child) in enumerate(element.children) + if (typeof(child) <: HTMLElement{:script}) && + haskey(child.attributes, "src") && + contains(child.attributes["src"], "data:application/javascript;base64,") + src_raw = replace( + child.attributes["src"], + "data:application/javascript;base64," => "", + ) + js_content = HTTP.base64decode(src_raw) + global js_counter += 1 + js_filename = "./js/js-$(js_counter).js" + js_write = "$dir/../_build/html/julia/js/js-$(js_counter).js" + write(js_write, js_content) + element.children[i] = HTMLElement{:script}( + Vector{HTMLNode}[], + element, + Dict{String,String}("src" => js_filename, "type" => "module"), + ) + end + if typeof(child) <: HTMLElement + separate_javascript!(child) + end + end + return element +end + +function separate_javascript(html_source::AbstractString) + if !isdir("$dir/../_build/html/") + mkdir("$dir/../_build/html/") + end + if !isdir("$dir/../_build/html/julia/") + mkdir("$dir/../_build/html/julia/") + end + if !isdir("$dir/../_build/html/julia/js/") + mkdir("$dir/../_build/html/julia/js/") + end + parsed_html = parsehtml(html_source) + div_element = separate_javascript!(parsed_html.root.children[2].children[1]) +end + +function Base.display(fig::FigureContainer{<:BackendType}) + display(fig.fig) +end + +function Base.display(fig::FigureContainer{WGLMakieType}) + html_str = string(DOM.div(fig.fig)) + html_source_new = separate_javascript(html_str) + display("text/html", string(html_source_new)) +end + +end From 47c2fd7f43253f727b6b022af26762be80dce1af Mon Sep 17 00:00:00 2001 From: Ian Hammond Date: Sun, 31 Mar 2024 15:24:04 -0400 Subject: [PATCH 7/8] initial nonlinear thermal --- docs/julia/heater_3d.jl | 46 ++++++++++++++--- src/.DS_Store | Bin 0 -> 6148 bytes src/Thermal/Thermal.jl | 110 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 148 insertions(+), 8 deletions(-) create mode 100644 src/.DS_Store diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl index dcdb8d0b..ea61cad2 100644 --- a/docs/julia/heater_3d.jl +++ b/docs/julia/heater_3d.jl @@ -244,7 +244,7 @@ electrical_conductivity = [ "metal3#e1" => 3.77e7, "metal3#e2" => 3.77e7, ] -thermal_conductivities = [ +thermal_conductivities_base = [ "core" => 90, "box" => 1.38, "clad" => 1.38, @@ -256,6 +256,18 @@ thermal_conductivities = [ "metal3#e1" => 400, "metal3#e2" => 400, ] +thermal_conductivities_linear= [ + "core" => 0.05*10, + "box" => 0.01*10, + "clad" => 0.01*10, + "heater" => 0.1*10, + "via2" => 0.1*10, + "metal2" => 0.1*10, + "via1" => 0.1*10, + "metal3" => 0.1*10, + "metal3#e1" => 0.1*10, + "metal3#e2" => 0.1*10, +] thermal_diffisitivities = [ "core" => 90 / 711 / 2330, "box" => 1.38 / 1000 / 2650, @@ -280,18 +292,22 @@ labels = get_face_labeling(model) tags = get_face_tag(labels, num_cell_dims(model)) τ = CellField(tags, Ω) +# Decide whether to use a nonlinear or linear thermal solve +nonlinear_flag = true + 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_conductivities_base = + Dict(get_tag_from_name(labels, u) => v for (u, v) in thermal_conductivities_base) +thermal_conductivities_linear = + Dict(get_tag_from_name(labels, u) => v for (u, v) in thermal_conductivities_linear) thermal_diffisitivities = Dict(get_tag_from_name(labels, u) => v for (u, v) in thermal_diffisitivities) ϵ_diffisitivities(tag) = thermal_diffisitivities[tag] + # %% [markdown] # The next step is to define the boundary conditions, this can be done simply via julia-dicts: @@ -313,7 +329,20 @@ println("Current: ", @sprintf("%.2f mA", current * 1e3)) println("Power: ", @sprintf("%.2f mW", power * 1e3)) # %% tags=["remove-stderr"] -T0 = calculate_temperature(ϵ_conductivities ∘ τ, power_density(p0), boundary_temperatures) +T0, ϵ_conductivities = if nonlinear_flag + ϵ_conductivities_base(tag) = thermal_conductivities_base[tag] + ϵ_conductivities_linear(tag) = thermal_conductivities_linear[tag] + ϵ_conductivities = u -> (ϵ_conductivities_base ∘ τ) + (ϵ_conductivities_linear ∘ τ) * u + Gridap.gradient(::typeof(ϵ_conductivities), du) = (ϵ_conductivities_linear ∘ τ) * du + + T0 = calculate_temperature(ϵ_conductivities, power_density(p0), boundary_temperatures; solver=NLSolver(method=:newton)) + T0, ϵ_conductivities +else + ϵ_conductivities = tag -> thermal_conductivities_base[tag] + + T0 = calculate_temperature(ϵ_conductivities ∘ τ, power_density(p0), boundary_temperatures; solver=BackslashSolver()) + T0, ϵ_conductivities ∘ τ +end # %% [markdown] # Yay, now we have both fields simulated! Let's get the average temperature of the silicon waveguide. @@ -405,7 +434,7 @@ ax = Axis( # %% tags=["remove-stderr", "hide-output"] Δt = 2e-7 t_end = 2e-4 -plot_every = 200 +plot_every = 250 timestamps = collect(Δt:Δt*plot_every:t_end) total_render_time = 5.0 fps = ceil(Int, length(timestamps) / total_render_time) @@ -413,13 +442,14 @@ html_sources = String[] lins = CairoMakie.Lines[] for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", 0, 1)] uₕₜ = calculate_temperature_transient( - ϵ_conductivities ∘ τ, + ϵ_conductivities, ϵ_diffisitivities ∘ τ, power_density(p0) * power_factor, boundary_temperatures, temperature(T0) * temperature_factor, Δt, t_end, + solver=ifelse(nonlinear_flag, NLSolver(method=:newton), LUSolver()), ) # Create animation diff --git a/src/.DS_Store b/src/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..84ebb0abc4c27b38370ddfe9bc915c6365164bf4 GIT binary patch literal 6148 zcmeHKy-veG47N)pg1VFm#(M+Qg*jAVU|=c_P)aL0B&0$@Fy1pU!pvjv3Ooa!&u&SS zR$@Y^vL)Yl@!5CImml8|5s#ndBcdS@Wl+K15t=O`1axfHU9>>>UHB*&^9M z(OYN08E^(R49Nc>Km~KdsF;2oXwne?*oQd@`f1Mu3}OK0hEWj~2x}-%L)nfPtl_W+ z^UDpRqJ|Tvl`)U4%yvRy+p%#-?!>vGx6Xhw&}QI3FNae9FL$5++d+Qi3^)UO#Q^t< zNioJPS#9myoYdL?okB&#FDh g(x, t) + return g + end + + tags_temperatures = [def_g(v) for (u, v) in temperatures_c] + if isempty(tags_temperatures) + tags_temperatures = x -> 0 + end + + model = get_active_model(get_triangulation(thermal_diffusitivity)) + V = TestFESpace( + model, + ReferenceFE(lagrangian, Float64, order); + conformity = :H1, + dirichlet_tags = tags, + ) + U = TransientTrialFESpace(V, tags_temperatures) + Ω = Triangulation(model) + dΩ = Measure(Ω, order) + + m₀(t, u, v) = ∫(thermal_conductivity(u) / thermal_diffusitivity * u * v)dΩ + a₀(t, u, v) = ∫(thermal_conductivity(u) * (∇(u) ⋅ ∇(v)))dΩ + + b = if typeof(power_density) <: Function + println("function!") + (t,v) -> ∫(power_density(t) * v)dΩ + else + (t,v) -> ∫(power_density * v)dΩ + end + + res(t, u, v) = m₀(t, u, v) + a₀(t, u, v) - b(t, v) + jac(t, u, du, v) = ∫( + ∇(v) ⋅ ( + (thermal_conductivity(u)) * ∇(du) + + ∇(thermal_conductivity)(du) * ∇(u) + ) + )*dΩ + jac_t(t,u,duₜ,v) = ∫( duₜ*v )dΩ + + op_C = TransientFEOperator(res,jac,jac_t,U,V) + + θ = 0.5 + ode_solver = ThetaMethod(solver, Δt, θ) + + @show "beginning interpolation" + u₀ = interpolate_everywhere(T0, U(0.0)) + @show "beginning solve" + uₕₜ = solve(ode_solver, op_C, u₀, 0, t_end) + @show "end solve" + + return uₕₜ +end + end From e7086b159cbc07a028ceb9931f163dca1a4df766 Mon Sep 17 00:00:00 2001 From: Ian Hammond Date: Sun, 31 Mar 2024 15:26:02 -0400 Subject: [PATCH 8/8] removed accidental print statements --- src/Thermal/Thermal.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/Thermal/Thermal.jl b/src/Thermal/Thermal.jl index 5a1b2bb5..3318dc4e 100644 --- a/src/Thermal/Thermal.jl +++ b/src/Thermal/Thermal.jl @@ -201,11 +201,8 @@ function calculate_temperature_transient( θ = 0.5 ode_solver = ThetaMethod(solver, Δt, θ) - @show "beginning interpolation" u₀ = interpolate_everywhere(T0, U(0.0)) - @show "beginning solve" uₕₜ = solve(ode_solver, op_C, u₀, 0, t_end) - @show "end solve" return uₕₜ end