@@ -3,7 +3,7 @@ using SSP: Kernel, Pad, Convolve, Project
33using . Kernel: conickernel
44using . Pad: FillPadding, BoundaryPadding, Inner, PaddingProblem, DefaultPaddingAlgorithm
55using . Convolve: DiscreteConvolutionProblem, FFTConvolution
6- using . Project: ProjectionProblem, SSP2
6+ using . Project: ProjectionProblem, SSP1_linear, SSP1, SSP2
77
88using Random
99using CairoMakie
@@ -29,8 +29,8 @@ kernel = conickernel(grid, radius)
2929
3030padprob = PaddingProblem (;
3131 data = design_vars,
32- # boundary = BoundaryPadding(size(kernel) .- 1, size(kernel) .- 1),
33- boundary = FillPadding (1.0 , size (kernel) .- 1 , size (kernel) .- 1 ),
32+ boundary = BoundaryPadding (size (kernel) .- 1 , size (kernel) .- 1 ),
33+ # boundary = FillPadding(1.0, size(kernel) .- 1, size(kernel) .- 1),
3434)
3535padalg = DefaultPaddingAlgorithm ()
3636padsolver = init (padprob, padalg)
@@ -57,8 +57,8 @@ filtered_design_vars = depadsol.value
5757
5858# projection points need not be the same as design variable grid
5959target_grid = (
60- range (- 1 , 1 , length= Nx * 2 ),
61- range (- 1 , 1 , length= Ny * 2 ),
60+ range (- 1 , 1 , length= Nx * 1 ),
61+ range (- 1 , 1 , length= Ny * 1 ),
6262)
6363target_points = vec (collect (Iterators. product (target_grid... )))
6464projprob = ProjectionProblem (;
@@ -68,6 +68,8 @@ projprob = ProjectionProblem(;
6868 beta = Inf ,
6969 eta = 0.5 ,
7070)
71+ # projalg = SSP1_linear()
72+ # projalg = SSP1()
7173projalg = SSP2 ()
7274projsolver = init (projprob, projalg)
7375projsol = solve! (projsolver)
8789end
8890
8991function fom (data, grid)
90- return sum (abs2, data) * prod (step, grid )
92+ return sum (abs2, data) / length (data )
9193end
9294obj = fom (projected_design_vars, grid)
9395
9496function adjoint_fom (adj_fom, data, grid)
9597 adjoint_fom! (similar (data), adj_fom, data, grid)
9698end
9799function adjoint_fom! (adj_data, adj_fom, data, grid)
98- adj_data .= adj_fom .* 2 .* data .* prod (step, grid)
100+ adj_data .= ( adj_fom / length (data)) .* 2 .* data
99101 return adj_data
100102end
101103
189191 projsol = solve! (projsolver)
190192
191193 fig = Figure ()
192- ax1 = Axis (fig[1 ,1 ]; title = " Objective history" , yscale= log10)
194+ ax1 = Axis (fig[1 ,1 ]; title = " Objective history" , yscale= log10, limits = ( nothing , ( 1e-16 , 1e1 )) )
193195 h1 = scatterlines! (ax1, evaluation_history)
194196
195197 ax2 = Axis (fig[1 ,2 ]; title = " Final SSP2 design" , aspect= DataAspect ())
0 commit comments