Hacktoberfest 2026: die Issues, die Maintainer für den Oktober markiert haben – offen und einsteigerfreundlich. Hacktoberfest-Issues durchsuchen

Issue in running Tutorial 18 - Topology optimization

Offen
  • Stimmt etwas nicht?
#180 0 Kommentare 0 Reaktionen 0 zugewiesene Personen Auf GitHub ansehen

Dieses Issue hat noch niemand übernommen.

Bewertung

Schwierigkeit
3/5
Geschätzter Aufwand
1-2 Tage
Anfängerfreundlichkeit
35/100
Issue-Typ
Bug
Klarheit
Größtenteils klar
Aktivitätsstatus
Veraltet
Tech-Stack
julia
Bereich
documentation

Rechercherichtung

Führen Sie zunächst Tutorial 18 - Topology optimization mit dem bereitgestellten Code aus und untersuchen Sie MatrixA, insbesondere a_design und ξd, bei denen der MethodError während der Assemblierung auftritt. Überprüfen Sie, dass das Tutorial über A_mat = MatrixA(pth; phys_params, fem_params) hinaus ausgeführt wird, ohne den ForwardDiff.Dual-to-Float64-Fehler.

Vom Indexierungsmodell aus dem Issue-Text verfasst.

Beschreibung

Hi,
I get the following error when running the Tutorial 18 - Topology optimization

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{Gridap.Fields.var"#89#90"{Λ{NamedTuple{(:k, :n_metal, :n_air, :μ, :R, :dpml, :LHp, :LHn), Tuple{Float64, ComplexF64, Int64, Int64, Float64, Int64, Tuple{Float64, Int64}, Tuple{Float64, Int64}}}}}, Float64}, Float64, 2})

The error appears to be at line A_mat = MatrixA(pth; phys_params, fem_params)

Please let me know how to fix it. Thanks!

using Gridap, Gridap.Geometry, Gridap.Fields, GridapGmsh
λ = 532      # Wavelength (nm)
L = 600      # Width of the numerical cell (excluding PML) (nm)
h1 = 600     # Height of the air region below the source (nm)
h2 = 200     # Height of the air region above the source (nm)
dpml = 300   # Thickness of the PML (nm)

n_metal = 0.054 + 3.429im # Silver refractive index at λ = 532 nm
n_air = 1    # Air refractive index
μ = 1        # Magnetic permeability
k = 2*π/λ

model = GmshDiscreteModel("RecCirGeometry.msh")

order = 1
reffe = ReferenceFE(lagrangian, Float64, order)
V = TestFESpace(model, reffe, dirichlet_tags = ["DirichletEdges", "DirichletNodes"], vector_type = Vector{ComplexF64})
U = V 

degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)

Γ_s = BoundaryTriangulation(model; tags = ["Source"]) # Source line
dΓ_s = Measure(Γ_s, degree)

Ω_d = Triangulation(model, tags="Design")
dΩ_d = Measure(Ω_d, degree)

Ω_c = Triangulation(model, tags="Center")
dΩ_c = Measure(Ω_c, degree)

p_reffe = ReferenceFE(lagrangian, Float64, 0)
Q = TestFESpace(Ω_d, p_reffe, vector_type = Vector{Float64})
P = Q
np = num_free_dofs(P)

pf_reffe = ReferenceFE(lagrangian, Float64, 1)
Qf = TestFESpace(Ω_d, pf_reffe, vector_type = Vector{Float64})
Pf = Qf
fem_params = (; V, U, Q, P, Qf, Pf, np, Ω, dΩ, dΩ_d, dΩ_c, dΓ_s)
R = 1e-10
LHp=(L/2, h1+h2)   # Start of PML for x,y > 0
LHn=(L/2, 0)       # Start of PML for x,y < 0
phys_params = (; k, n_metal, n_air, μ, R, dpml, LHp, LHn)

function s_PML(x; phys_params)
    σ = -3 / 4 * log(phys_params.R) / phys_params.dpml / phys_params.n_air
    xf = Tuple(x)
    u = @. ifelse(xf > 0 , xf - phys_params.LHp, - xf - phys_params.LHn)
    return @. ifelse(u > 0,  1 + (1im * σ / phys_params.k) * (u / phys_params.dpml)^2, $(1.0+0im))
end

function ds_PML(x; phys_params)
    σ = -3 / 4 * log(phys_params.R) / phys_params.dpml / phys_params.n_air
    xf = Tuple(x)
    u = @. ifelse(xf > 0 , xf - phys_params.LHp, - xf - phys_params.LHn)
    ds = @. ifelse(u > 0, (2im * σ / phys_params.k) * (1 / phys_params.dpml)^2 * u, $(0.0+0im))
    return ds.*sign.(xf)
end

struct Λ{PT} <: Function
    phys_params::PT
end

function (Λf::Λ)(x)
    s_x,s_y = s_PML(x; Λf.phys_params)
    return VectorValue(1/s_x, 1/s_y)
end

r = 5/sqrt(3)               # Filter radius
β = 32.0                    # β∈[1,∞], threshold sharpness
η = 0.5                     # η∈[0,1], threshold center

a_f(r, u, v) = r^2 * (∇(v) ⋅ ∇(u))

function Filter(p0; r, fem_params)
    ph = FEFunction(fem_params.P, p0)
    op = AffineFEOperator(fem_params.Pf, fem_params.Qf) do u, v
        ∫(a_f(r, u, v))fem_params.dΩ_d + ∫(v * u)fem_params.dΩ_d, ∫(v * ph)fem_params.dΩ_d
      end
    pfh = solve(op)
    return get_free_dof_values(pfh)
end

function Threshold(pfh; β, η)
    return ((tanh(β * η) + tanh(β * (pfh - η))) / (tanh(β * η) + tanh(β * (1.0 - η))))
end

using LinearAlgebra
ξd(p, n_air, n_metal)= 1 / (n_air + (n_metal - n_air) * p)^2 - 1 / n_air^2 # in the design region

a_base(u, v; phys_params) = (1 / phys_params.n_air^2) * ((∇ .* (Λ(phys_params) * v)) ⊙ (Λ(phys_params) .* ∇(u))) - (phys_params.k^2 * phys_params.μ * (v * u))

a_design(u, v, pth; phys_params) = ((p -> ξd(p, phys_params.n_air, phys_params.n_metal)) ∘ pth) * (∇(v) ⊙ ∇(u))

function MatrixA(pth; phys_params, fem_params)
    A_mat = assemble_matrix(fem_params.U, fem_params.V) do u, v
        ∫(a_base(u, v; phys_params))fem_params.dΩ + ∫(a_design(u, v, pth; phys_params))fem_params.dΩ_d
    end
    return lu(A_mat)
end

p0 = zeros(fem_params.np)  # Here we make p=0 everywhere just for illustration purpose
pf_vec = Filter(p0;r, fem_params)
pfh = FEFunction(fem_params.Pf, pf_vec)
pth = (pf -> Threshold(pf; β, η)) ∘ pfh
A_mat = MatrixA(pth; phys_params, fem_params)
b_vec = assemble_vector(v->(∫(v)fem_params.dΓ_s), fem_params.V)
u_vec = A_mat \ b_vec
uh = FEFunction(fem_params.U, u_vec)
Vorherrschende Sprache
Julia
Sterne
145
Forks
50
Keine gemergten PRs in 30 T.

Entwicklungsumgebung

Erste Schritte

  1. Lesen Sie das ganze Issue und danach den Beitragsleitfaden des Projekts.
  2. Schreiben Sie ins Issue, dass Sie es übernehmen — das erspart doppelte Arbeit.
  3. Forken Sie das Repository und arbeiten Sie in einem Branch.
  4. Öffnen Sie einen Pull Request, der die Issue-Nummer nennt.

Mehr aus gridap/Tutorials

Alle Issues in gridap/Tutorials

Ähnliche Issues

Weitere Issues zu Julia

Neue Issues direkt in Ihr Postfach

Eine kurze Übersicht über anfängerfreundliche GitHub-Issues.