Skip to content

Commit

Permalink
Avoid adding empty particles blowing up ExaFMM; make sure unit tests …
Browse files Browse the repository at this point in the history 
…are passing
  • Loading branch information
@EdoAlvarezR
EdoAlvarezR committed May 13, 2021
1 parent 9077eae commit d302a74
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 48 deletions.
103 changes: 56 additions & 47 deletions examples/tandemheavingwing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,41 +226,51 @@ function tandemheavingwing(; # TEST OPTIONS
prev_wing = nothing
prev_system = nothing

fig1, fig2 = nothing, nothing
axs1, axs2 = nothing, nothing


function monitor(sim, PFIELD, T, DT; figname="monitor_$(save_path)",
nsteps_plot=1)

aux = PFIELD.nt/nsteps
clr = (1-aux, 0, aux)

if PFIELD.nt==0 && disp_plot
figure(figname, figsize=[7*2, 5*2]*figsize_factor)
subplot(221)
xlim([0,1])
xlabel(L"$\frac{2y}{b}$")
ylabel(L"$\frac{Cl}{CL}$")
title("Spanwise lift distribution")

subplot(222)
xlim([0,1])
xlabel(L"$\frac{2y}{b}$")
ylabel(L"$\frac{Cd}{CD}$")
title("Spanwise drag distribution")

subplot(223)
xlabel("Simulation time (s)")
ylabel(L"Lift Coefficient $C_L$")

subplot(224)
xlabel("Simulation time (s)")
ylabel(L"Drag Coefficient $C_D$")

figure(figname*"_2", figsize=[7*2, 5*1]*figsize_factor)
subplot(121)
xlabel(L"$\frac{2y}{b}$")
ylabel(L"Circulation $\Gamma$")
subplot(122)
xlabel(L"$\frac{2y}{b}$")
ylabel(L"Effective velocity $V_\infty$")
fig1 = figure(figname, figsize=[7*2, 5*2]*figsize_factor)
axs1 = fig1.subplots(2, 2)
axs1 = [axs1[1], axs1[3], axs1[2], axs1[4]]

ax = axs1[1]
ax.set_xlim([0,1])
ax.set_xlabel(L"$\frac{2y}{b}$")
ax.set_ylabel(L"$\frac{Cl}{CL}$")
ax.set_title("Spanwise lift distribution")

ax = axs1[2]
ax.set_xlim([0,1])
ax.set_xlabel(L"$\frac{2y}{b}$")
ax.set_ylabel(L"$\frac{Cd}{CD}$")
ax.set_title("Spanwise drag distribution")

ax = axs1[3]
ax.set_xlabel("Simulation time (s)")
ax.set_ylabel(L"Lift Coefficient $C_L$")

ax = axs1[4]
ax.set_xlabel("Simulation time (s)")
ax.set_ylabel(L"Drag Coefficient $C_D$")

fig2 = figure(figname*"_2", figsize=[7*2, 5*1]*figsize_factor)
axs2 = fig2.subplots(1, 2)

ax = axs2[1]
ax.set_xlabel(L"$\frac{2y}{b}$")
ax.set_ylabel(L"Circulation $\Gamma$")

ax = axs2[2]
ax.set_xlabel(L"$\frac{2y}{b}$")
ax.set_ylabel(L"Effective velocity $V_\infty$")

end

Expand Down Expand Up @@ -296,30 +306,29 @@ function tandemheavingwing(; # TEST OPTIONS
vlm._addsolution(wing, "Cl/CL", ClCL)
vlm._addsolution(wing, "Cd/CD", CdCD)

subplot(221)
plot(web_2yb, web_ClCL, "ok", label="Weber's experimental data")
plot(y2b, ClCL, "-", label="FLOWVLM", alpha=0.5, color=clr)
ax = axs1[1]
ax.plot(web_2yb, web_ClCL, "ok", label="Weber's experimental data")
ax.plot(y2b, ClCL, "-", label="FLOWVLM", alpha=0.5, color=clr)

subplot(222)
plot(web_2yb, web_CdCD, "ok", label="Weber's experimental data")
plot(y2b, CdCD, "-", label="FLOWVLM", alpha=0.5, color=clr)
ax = axs1[2]
ax.plot(web_2yb, web_CdCD, "ok", label="Weber's experimental data")
ax.plot(y2b, CdCD, "-", label="FLOWVLM", alpha=0.5, color=clr)

subplot(223)
plot([0, T], web_CL*ones(2), ":k", label="Weber's experimental data")
plot([T], [CLwing], "o", label="FLOWVLM", alpha=0.5, color=clr)
ax = axs1[3]
ax.plot([0, T], web_CL*ones(2), ":k", label="Weber's experimental data")
ax.plot([T], [CLwing], "o", label="FLOWVLM", alpha=0.5, color=clr)

subplot(224)
plot([0, T], web_CD*ones(2), ":k", label="Weber's experimental data")
plot([T], [CDwing], "o", label="FLOWVLM", alpha=0.5, color=clr)
ax = axs1[4]
ax.plot([0, T], web_CD*ones(2), ":k", label="Weber's experimental data")
ax.plot([T], [CDwing], "o", label="FLOWVLM", alpha=0.5, color=clr)

figure(figname*"_2")
subplot(121)
plot(y2b, wing.sol["Gamma"], "-", label="FLOWVLM", alpha=0.5, color=clr)
ax = axs2[1]
ax.plot(y2b, wing.sol["Gamma"], "-", label="FLOWVLM", alpha=0.5, color=clr)
if wake_coupled && PFIELD.nt!=0
subplot(122)
plot(y2b, norm.(wing.sol["Vkin"])/magVinf, "-", label="FLOWVLM", alpha=0.5, color=[clr[1], 1, clr[3]])
if VehicleType!=uns.QVLMVehicle; plot(y2b, norm.(wing.sol["Vvpm"]), "-", label="FLOWVLM", alpha=0.5, color=clr); end;
plot(y2b, [norm(Vinf(vlm.getControlPoint(wing, i), T)) for i in 1:vlm.get_m(wing)],
ax = axs2[2]
ax.plot(y2b, norm.(wing.sol["Vkin"])/magVinf, "-", label="FLOWVLM", alpha=0.5, color=[clr[1], 1, clr[3]])
if VehicleType!=uns.QVLMVehicle; ax.plot(y2b, norm.(wing.sol["Vvpm"]), "-", label="FLOWVLM", alpha=0.5, color=clr); end;
ax.plot(y2b, [norm(Vinf(vlm.getControlPoint(wing, i), T)) for i in 1:vlm.get_m(wing)],
"-k", label="FLOWVLM", alpha=0.5)
end

Expand Down
8 changes: 7 additions & 1 deletion src/FLOWUnsteady_simulation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,12 @@ function add_particle(pfield::vpm.ParticleField, X::Array{Float64, 1},

Gamma = gamma*(V*dt)*infD # Vectorial circulation

# Avoid adding empty particles to the computational domain, or ExaFMM would
# blow up
if sqrt(Gamma[1]^2 + Gamma[2]^2 + Gamma[3]^2) <= 5*eps(vpm.RealFMM)
Gamma = 5*eps(vpm.RealFMM)*ones(3)
end

# Decreases p_per_step for slowly moving parts of blade
# aux = min((sigma/p_per_step)/overwrite_sigma, 1)
# pps = max(1, min(p_per_step, floor(Int, 1/(1-(aux-1e-14))) ))
Expand All @@ -229,7 +235,7 @@ function add_particle(pfield::vpm.ParticleField, X::Array{Float64, 1},
dX = l/pps
for i in 1:pps
vpm.add_particle(pfield, X + i*dX - dX/2, Gamma/pps, sigmap;
vol=vol/pps, circulation=abs(gamma))
vol=vol/pps, circulation=max(abs(gamma), 5*eps(vpm.RealFMM)))
end
end

Expand Down

0 comments on commit d302a74

Please sign in to comment.