Open masmoudiept opened 3 months ago
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)
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)
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!