Hacktoberfest 2026 : les issues que les mainteneurs ont marquées pour octobre, ouvertes et accessibles aux débutants. Parcourir les issues Hacktoberfest

Issue in running Tutorial 18 - Topology optimization

Ouverte
  • Quelque chose ne va pas ?
#180 0 commentaires 0 réactions 0 personnes assignées Voir sur GitHub

Personne n'a encore pris cette issue.

Évaluation

Difficulté
3/5
Temps estimé
1-2 jours
Accessibilité débutants
35/100
Type d'issue
Bug
Clarté
Plutôt claire
Activité
À l'abandon
Stack technique
julia
Domaine
documentation

Piste de recherche

Commencez par exécuter Tutorial 18 - Topology optimization avec le code fourni et inspectez MatrixA, en particulier a_design et ξd, où le MethodError apparaît pendant l’assemblage. Vérifiez que le tutoriel s’exécute au-delà de A_mat = MatrixA(pth; phys_params, fem_params) sans l’erreur ForwardDiff.Dual-to-Float64.

Rédigé par le modèle d'indexation à partir du texte de l'issue.

Description

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)
Langage dominant
Julia
Étoiles
145
Forks
50
Aucune PR mergée en 30 j

Préparer son environnement

Par où commencer

  1. Lisez l'issue en entier, puis le guide de contribution du projet.
  2. Signalez en commentaire que vous la prenez — cela évite que deux personnes fassent le même travail.
  3. Forkez le dépôt et travaillez sur une branche.
  4. Ouvrez une pull request qui référence le numéro de l'issue.

Autres issues de gridap/Tutorials

Toutes les issues de gridap/Tutorials

Issues similaires

Plus d'issues Julia

Recevez les nouvelles issues par e-mail

Un résumé court des issues GitHub adaptées aux débutants.