Skip to content
Snippets Groups Projects
Commit e9dd08534a05 authored by Pierre Augier's avatar Pierre Augier
Browse files

py/save_forcing_vs_time.py

parent d8dfddf1146e
No related branches found
No related tags found
1 merge request!12Topic/default/update strat
......@@ -35,6 +35,7 @@
print(f"Calling {source} to make {name}")
run_path(str(here / source))
make_fig("forcing_vs_time")
make_table("better_simuls")
......@@ -38,21 +39,21 @@
make_table("better_simuls")
make_fig(f"figure2")
make_fig(f"figure3")
make_fig(f"figure5")
make_fig(f"figure6")
make_fig(f"figure7")
make_fig(f"figure8")
make_fig(f"figure9")
make_fig(f"figure10")
make_fig(f"figure11")
make_fig(f"figure12")
make_fig(f"figure13")
make_fig(f"figure14")
make_fig(f"figure15")
make_fig(f"figure16")
make_fig(f"figure17")
make_fig("figure2")
make_fig("figure3")
make_fig("figure5")
make_fig("figure6")
make_fig("figure7")
make_fig("figure8")
make_fig("figure9")
make_fig("figure10")
make_fig("figure11")
make_fig("figure12")
make_fig("figure13")
make_fig("figure14")
make_fig("figure15")
make_fig("figure16")
make_fig("figure17")
"""
for n in (1, 17):
......
import numpy as np
import matplotlib.pyplot as plt
import scipy
from fluiddyn.io.redirect_stdout import stdout_redirected
from fluidsim.solvers.ns3d.strat.solver import Simul
from util import save_fig
params = Simul.create_default_params()
params.output.HAS_TO_SAVE = False
params.output.sub_directory = "tmp4fig"
params.forcing.tcrandom.time_correlation = 1.0
params.oper.nx = nx = 12
params.oper.ny = params.oper.nz = nx
params.init_fields.type = "noise"
params.init_fields.noise.length = params.oper.Lz / 2
params.init_fields.noise.velo_max = 1.0
params.forcing.nkmin_forcing = 2
params.forcing.nkmax_forcing = nx // 4
params.forcing.enable = True
params.forcing.type = "tcrandom"
params.forcing.normalized.type = "2nd_degree_eq"
# this actually produces a cleaner forcing...
# params.forcing.normalized.which_root = "positive"
params.forcing.key_forced = "vx_fft"
with stdout_redirected():
sim = Simul(params)
dt = 0.01
t_end = 50
freq_sample = 1 / dt
nt = int(t_end / dt)
fx_vs_time = np.empty([nx**3, nt])
times = []
# effect of the forcing on the state
epsilon = 2e-3
dissipation_coef = 1 - 1e-3
for idx in range(nt):
time = idx * dt
times.append(time)
sim.time_stepping.t = time
sim.forcing.compute()
forcing_spect = sim.forcing.get_forcing()
sim.state.state_spect = (
dissipation_coef * sim.state.state_spect + epsilon * forcing_spect
)
fx_fft = forcing_spect.get_var("vx_fft")
fx = sim.oper.ifft(fx_fft).flatten()
fx_vs_time[:, idx] = fx
times = np.array(times)
freqs, periodograms = scipy.signal.periodogram(
fx_vs_time, freq_sample, scaling="density"
)
periodogram = periodograms.mean(0)
periodogram[0] = np.nan
fig, (ax0, ax1) = plt.subplots(2)
indices = list(range(0, 80, 10))
tmin_plot = 20
tmax_plot = tmin_plot + 10
# tmax_plot = max(times)
# tmin_plot = tmax_plot - 10
index_tmax = abs(times - tmax_plot).argmin()
index_tmin = abs(times - tmin_plot).argmin()
times = times[index_tmin : index_tmax + 1]
fx_vs_time = fx_vs_time[indices, index_tmin : index_tmax + 1]
for ix in range(fx_vs_time.shape[0]):
ax0.plot(times, fx_vs_time[ix])
ax0.set_xlabel("time / correlation time")
ax0.set_title("forcing for few spatial points")
ax0.set_xlim((tmin_plot, tmax_plot))
low, high = ax0.get_ylim()
bound = max(abs(low), abs(high))
ax0.set_ylim(-bound, bound)
ax1.loglog(freqs, periodogram)
ax1.set_xlabel("frequency / correlation frequency")
ax1.set_ylabel("PDS forcing")
ax1.set_xlim(right=10)
ax1.set_ylim(bottom=1e-4)
fig.tight_layout()
save_fig(fig, "forcing_vs_time.png")
......@@ -15,9 +15,5 @@
path_base = os.environ["STRAT_TURB_POLO2022"]
path_base_proj = os.environ["STRAT_TURB_POLO_PROJ2022"]
path_base_ratio_one = os.environ["STRAT_WAVES2022"]
paths_all = sorted(Path(path_base).glob("simul_folders/ns3d*"))
paths_all_proj = sorted(Path(path_base_proj).glob("simul_folders/ns3d*"))
......@@ -22,8 +18,14 @@
paths_all = sorted(Path(path_base).glob("simul_folders/ns3d*"))
paths_all_proj = sorted(Path(path_base_proj).glob("simul_folders/ns3d*"))
paths_all_ratio_one = sorted(
Path(path_base_ratio_one).glob("simul_folders/ns3d*")
)
try:
path_base_ratio_one = os.environ["STRAT_WAVES2022"]
except KeyError:
print("warning: STRAT_WAVES2022 environment variable not set")
else:
paths_all_ratio_one = sorted(
Path(path_base_ratio_one).glob("simul_folders/ns3d*")
)
here = Path(__file__).absolute().parent
tmp_dir = here.parent / "tmp"
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment