Skip to content

Commit

Permalink
add forward inverse cond plots
Browse files Browse the repository at this point in the history
  • Loading branch information
@arturgower
arturgower committed Nov 5, 2024
1 parent 194263d commit 5a09bef
Show file tree
Hide file tree
Showing 5 changed files with 135 additions and 40 deletions.
106 changes: 86 additions & 20 deletions docs/examples/condition_numbers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,17 +8,17 @@ medium = Elastic(2; ρ = 7000.0, cp = 5000.0 - 0.0im, cs = 3500.0 - 0.0im)

inner_radius = 1.0
bearing = RollerBearing(medium = medium,
inner_radius = inner_radius, outer_radius = inner_radius + 0.1,
angular_speed = Ω,
rollers_inside = true,
number_of_rollers = Z,
roller_radius = inner_radius / (typeof(inner_radius)(1.0) + 2.5*Z / (2pi))
inner_radius = inner_radius,
outer_radius = inner_radius + 0.3,
)

krs = LinRange(0,50,400);
krs = LinRange(0,25,300);
# krs = LinRange(0,5,240);
ωs = real(medium.cp) .* krs

ns = 0:50
ns = 0:40
ns = 0:20

# Types of boundary conditions for the forward and inverse problem
bc1_forward = TractionBoundary(inner=true)
Expand All @@ -42,7 +42,8 @@ for n in ns, ω in ωs]
using LaTeXStrings

# pyplot(size = (300,250))
gr(size = (300,250))
gr(size = (250,250))
gr(size = (310,250))
# Plots.PyPlot.rc("text", usetex="true")

heatmap(krs, ns, conds, clims = (0,25),
Expand All @@ -52,12 +53,24 @@ heatmap(krs, ns, conds, clims = (0,25),

# savefig("docs/images/condition-forward-steel-r1-1.0-r2-1p5.pdf")

dkr = 8.0;
dn = 4

dkrs = 2.5:dkr:(2dkr+2.5) |> collect
dns = 2.0:dn:(dn+2) |> collect

dns2 = [n for kr in dkrs, n in dns]
dkrs2 = [kr for kr in dkrs, n in dns]

scatter!([kr for n in dns, kr in dkrs][:], [n for n in dns, kr in dkrs][:], lab = "")

# plot!(10:40, (10:40) .* 0.0 .+ 20.)
# plot!(xlims = (30,31), ylims = (17,21), clims = (0,7525))
plot!(xlims = (0,20), ylims = (0,20), aspect_ratio = 1.0)
# savefig("docs/images/condition-scat-forward-steel-r1-1.0-r2-1p3.pdf")


bc1_inverse = DisplacementBoundary(inner=true)
bc2_inverse = TractionBoundary(inner=true)
bc1_inverse = DisplacementBoundary(outer=true)
bc2_inverse = TractionBoundary(outer=true)

conds = [
begin
Expand All @@ -81,24 +94,77 @@ heatmap(krs, ns, conds, clims = (0,25),
# savefig("docs/images/condition-inverse-steel.pdf")


bc1_inverse = DisplacementBoundary(inner=true)
bc2_inverse = DisplacementBoundary(outer=true)
## Check forward and inverse solve the same problem

conds = [
krs = LinRange(0,50,500);
ωs = real(medium.cp) .* krs
ns = 0:50

n = -3;
ω = ωs[10]

δ = 1e-8
method = ModalMethod(regularisation_parameter = δ,
only_stable_modes = false
)

errors = [
begin
try
nondim_bearing = nondimensionalise(bearing, ω)
A = boundarycondition_system(ω, nondim_bearing, bc1_inverse, bc2_inverse, n)
SM = diagm([(4.0) / sum(abs.(A[:,j])) for j in 1:size(A,2)])
A = A * SM
cond(A)
f0 = [1.0 1.0];
amp = 0.02 .* mean(abs.(f0));
error0 = amp .* ((rand(1,2) .- 0.5) + (rand(1,2) .- 0.5) .* im)
error1 = amp .* ((rand(1,2) .- 0.5) + (rand(1,2) .- 0.5) .* im)
error1 = error1 .* 0.0

bd1 = BoundaryData(bc1_forward; coefficients = f0 + error0, modes = [n])
bd2 = BoundaryData(bc2_forward; coefficients = [0.0 0.0] + error1, modes = [n])

sim = BearingSimulation(ω, bearing, bd1, bd2;
method = method,
nondimensionalise = true
)
wave = ElasticWave(sim);

f1 = field_modes(wave, bearing.outer_radius, bc1_inverse.fieldtype)
f2 = field_modes(wave, bearing.outer_radius, bc2_inverse.fieldtype)

amp = 0.02 .* mean(abs.(f1));
f1 = f1 + amp .* ((rand(1,2) .- 0.5) + (rand(1,2) .- 0.5) .* im)

amp = 0.02 .* mean(abs.(f2));
f2 = f2 + amp .* ((rand(1,2) .- 0.5) + (rand(1,2) .- 0.5) .* im)

bd1 = BoundaryData(bc1_inverse; coefficients = f1, modes = [n])
bd2 = BoundaryData(bc2_inverse; coefficients = f2, modes = [n])

sim = BearingSimulation(ω, bearing, bd1, bd2;
method = method,
nondimensionalise = true
)
wave = ElasticWave(sim);

f1 = field_modes(wave, bearing.inner_radius, bc1_forward.fieldtype)
f2 = field_modes(wave, bearing.outer_radius, bc2_forward.fieldtype)

norm(f1 - f0) / sqrt(2)
# norm(f2 - [0.0 0.0]) / sqrt(2)
catch
Inf
end
end
for n in ns, ω in ωs]

heatmap(krs, ns, conds, clims = (0,25),
c = cgrad(:navia, rev = true)
# gr(size = (360,200))
gr(size = (310,250))

colorbar_ticks=(0:0.1:0.5, string.(round.(Int, (0:0.1:0.5) .* 100), "%"))

heatmap(krs, ns, errors, clims = (0,0.2), xlims = (0,50), ylims = (0,50)
, colorbar_ticks = colorbar_ticks
, aspect_ratio = 1.0
, c = cgrad(:navia, rev = true)
, xlab = L"k_p r_1"
, ylab = L"n")

# savefig("docs/images/forward-inverse-cond-r1-1p0-r1-1p1.pdf")
69 changes: 49 additions & 20 deletions docs/examples/plot-modes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,8 @@ medium = Elastic(2; ρ = 7000.0, cp = 5000.0 - 0.0im, cs = 3500.0 - 0.0im)

inner_radius = 1.0
bearing = RollerBearing(medium = medium,
inner_radius = inner_radius, outer_radius = inner_radius * 1.1,
angular_speed = Ω,
rollers_inside = true,
number_of_rollers = Z,
roller_radius = inner_radius / (typeof(inner_radius)(1.0) + 2.5*Z / (2pi))
inner_radius = inner_radius, outer_radius = inner_radius + 0.3,
number_of_rollers = 0,
)

krs = LinRange(15,40,24400);
Expand Down Expand Up @@ -61,34 +58,66 @@ bd2 = BoundaryData(TractionBoundary(outer=true); coefficients = 1e-4 .* [1.0 1.
fs = [f[1] for f in result_traction.field];
result_pressure = FrequencySimulationResult(fs, result_traction.x, [ω]);

plot(result_pressure, ω;
plot(result_pressure, result_pressure.ω[1];
seriestype=:heatmap
, field_apply = f -> real(f[1])
# , leg = false,
);
plot!(
xlims = (-bearing.inner_radius * 0.4, bearing.inner_radius * 0.4)
, ylims = (bearing.inner_radius * 0.8, bearing.outer_radius * 1.1)
,frame = :none, title="", xguide ="",yguide =""
# xlims = (-bearing.inner_radius * 0.4, bearing.inner_radius * 0.4),
# ylims = (bearing.inner_radius * 0.8, bearing.outer_radius * 1.1),
frame = :none, title="", xguide ="",yguide =""
)

dkr = 8.0;
dn = 4

dkrs = 2.5:dkr:(2dkr+2.5) |> collect
dns = 2.0:dn:(dn+2) |> collect

dns2 = reverse(Int[n for kr in dkrs, n in dns])
dkrs2 = [kr for kr in dkrs, n in dns]


results = map(eachindex(dkrs2)) do i

ω = real(medium.cp) * dkrs2[i];

bd1 = BoundaryData(TractionBoundary(inner=true); coefficients = [0.0 0.0], modes = [dns2[i]])
bd2 = BoundaryData(TractionBoundary(outer=true); coefficients = [1.0 1.0], modes = [dns2[i]])

sims = map(eachindex(ωs)) do i
# the option only_stable_modes = false means the method will try to solve for modes which are ill posed
method = ModalMethod(regularisation_parameter = δs[i],
tol = tol,
method = ModalMethod(
tol = tol,
only_stable_modes = false
)
BearingSimulation(ωs[i], bearing, bd1, bd2;
sim = BearingSimulation(ω, bearing, bd1, bd2;
method = method,
nondimensionalise = true
)
wave = ElasticWave(sim)

result_traction = field(wave, bearing, TractionType(); res = 250);
fs = [f[1] for f in result_traction.field];
result_pressure = FrequencySimulationResult(fs, result_traction.x, [ω]);
end

waves = [ElasticWave(s) for s in sims];
gr(size = (1200,800))

results = map(eachindex(ωs)) do i
result_traction = field(inverse_waves[i], bearing, TractionType(); res = 250);
fs = [f[1] for f in result_traction.field];
result_pressure = FrequencySimulationResult(fs, result_traction.x, [ωs[i]]);
end
ps = map(results) do r
plot(r, r.ω[1];
seriestype=:heatmap
, field_apply = f -> real(f[1])
, leg = false,
);
plot!(
# xlims = (-bearing.inner_radius * 0.4, bearing.inner_radius * 0.4),
# ylims = (bearing.inner_radius * 0.8, bearing.outer_radius * 1.1),
frame = :none, title="", xguide ="",yguide =""
)
plot!(bearing)
end;

gr(size = (1200,1200))
plot(ps..., layout = (2, 3))

# savefig("docs/images/modes-r1-1.0-r2-1p3.pdf")
Binary file added docs/images/condition-scat-forward-steel-r1-1.0-r2-1p3.pdf
View file
Binary file not shown.
Binary file added docs/images/forward-inverse-cond-r1-1p0-r1-1p1.pdf
View file
Binary file not shown.
Binary file added docs/images/modes-r1-1.0-r2-1p3.pdf
View file
Binary file not shown.

0 comments on commit 5a09bef

Please sign in to comment.