From e7a5915c66b33b6b64908e9c4cb5b2e73848ae91 Mon Sep 17 00:00:00 2001 From: paugier Date: Wed, 23 Sep 2020 00:31:53 +0200 Subject: [PATCH] Black for files in doc, bench, scripts --- Makefile | 2 +- bench/dedalus/ns2d.py | 42 ++-- bench/dedalus/ns2d_rot.py | 45 ++-- bench/dedalus/ns2d_rot_faster.py | 47 +++-- bench/dedalus/ns2dstrat.py | 67 +++--- bench/dedalus/ns2dstrat_rot.py | 76 +++---- bench/old/compare_fluidfft.py | 67 +++--- bench/old/plot_bench.py | 24 +-- bench/old/simul_bench.py | 6 +- bench/old/submit_benchs.py | 10 +- bench/profile_time_stepping/init.py | 9 +- bench/profiling/simul_profile2d.py | 4 +- bench/profiling/simul_profile_ns2d.py | 2 +- bench/profiling/simul_profile_ns2dstrat.py | 2 +- bench/profiling/simul_profile_ns3d.py | 3 +- bench/profiling/simul_profile_plate2d.py | 2 +- bench/profiling/simul_profile_sw1l.py | 2 +- bench/profiling/util_bench.py | 186 ++++++++++------- bench/spectraldns/bench_2d.py | 24 +-- bench/spectraldns/bench_3d.py | 81 ++++--- bench/spectraldns/profile_2d.py | 37 ++-- bench/spectraldns/profile_3d.py | 90 ++++---- bench/submit_benchmarks_legi.py | 29 +-- bench/util.py | 197 ++++++++++++------ doc/example_from_numpy_doc.py | 5 +- .../simul_idempotent.py | 12 +- .../simul_ns3dstrat_waves.py | 6 +- .../simul_ns3dstrat_waves_forcinginscript.py | 4 +- .../internal_waves_coriolis/submit_simul.py | 8 +- doc/examples/launch_simul.py | 10 +- doc/examples/milestone/run_simul.py | 6 +- doc/examples/simul_ad1d.py | 15 +- doc/examples/simul_ns2d.py | 12 +- doc/examples/simul_ns2d_plot.py | 19 +- doc/examples/simul_ns2d_plotphys.py | 23 +- doc/examples/simul_ns2d_restart.py | 1 - .../simul_ns2dbouss_initfields_in_script.py | 25 +-- doc/examples/simul_ns2dbouss_staticlayer.py | 47 +++-- .../simul_ns2dstrat_forcing_in_script.py | 42 ++-- doc/examples/simul_ns3dbouss_staticlayer.py | 78 +++---- .../simul_ns3dstrat_check_energy_conserv.py | 14 +- .../simul_ns3dstrat_forcing_random.py | 3 +- doc/examples/time_stepping/3d/run_simul.py | 20 +- .../bugs/no_energy_conservation_ns2dstrat.py | 8 +- scripts/bugs/wrong_output_it0_initfield.py | 12 +- scripts/launch/simul_test_conserveE.py | 46 ++-- scripts/launch/sw1l_force_vortex_grid.py | 43 ++-- scripts/launch/sw1l_no_flow.py | 51 +++-- scripts/ns2d.strat/coeff_diss.py | 195 ++++++++++------- scripts/ns2d.strat/compute_anisotropy.py | 15 +- scripts/ns2d.strat/compute_flow_features.py | 33 +-- scripts/ns2d.strat/compute_length_scales.py | 12 +- .../ns2d.strat/compute_ratio_dissipation.py | 23 +- scripts/ns2d.strat/compute_reynolds_froude.py | 29 +-- scripts/ns2d.strat/compute_spectrum_kykx.py | 89 ++++---- scripts/ns2d.strat/compute_vertical_length.py | 12 +- scripts/ns2d.strat/flow_features.py | 44 ++-- scripts/ns2d.strat/make_table_parameters.py | 43 ++-- .../ns2d.strat/make_video_spectrum_kykx.py | 120 ++++++----- scripts/ns2d.strat/ns2dstrat_lmode.py | 53 +++-- .../plot_buoyancy_reynolds_anisotropy.py | 98 ++++++--- .../plot_buoyancy_reynolds_dissipation.py | 82 +++++--- .../plot_buoyancy_reynolds_froude.py | 80 +++++-- scripts/ns2d.strat/plot_froude_anisotropy.py | 96 ++++++--- scripts/ns2d.strat/plot_froude_dissipation.py | 82 +++++--- scripts/ns2d.strat/plot_length_scales.py | 92 +++++--- scripts/ns2d.strat/plot_reynolds_froude.py | 93 +++++++-- scripts/ns2d.strat/plot_vertical_scale.py | 99 ++++++--- .../ns2d.strat/plot_vertical_scale_godoy.py | 101 ++++++--- scripts/ns2d.strat/simul_from_state.py | 2 +- scripts/ns2d.strat/simul_init_linear_mode.py | 8 +- scripts/ns2d.strat/submit_legi.py | 10 +- scripts/ns2d.strat/submit_simul_from_state.py | 16 +- .../check_freq_spectra_with_linear_mode.py | 57 +++-- .../dispersion_relationship_SW1l.py | 40 ++-- scripts/plot_results/fig_spectra_forcingq.py | 114 ++++------ .../fig_spectra_forcingq_diff_c.py | 100 ++++----- scripts/plot_results/plot_interm.py | 171 ++++++++------- scripts/plot_results/plot_many_runs.py | 46 ++-- scripts/plot_results/plot_many_things.py | 20 +- scripts/plot_results/plot_profil_phys.py | 107 +++++----- scripts/plot_results/plot_spatial_means.py | 165 ++++++++------- scripts/plot_results/plot_time_means.py | 56 +++-- scripts/plot_results/superp_2Dspectra.py | 21 +- scripts/sphere.sw1l/galewsky2004.py | 22 +- scripts/util/modif_resol_all_dir.py | 12 +- scripts/util/modif_resolution.py | 19 +- scripts/util/resume_from_path.py | 7 +- 88 files changed, 2340 insertions(+), 1708 deletions(-) diff --git a/Makefile b/Makefile index c9b15a0f..a8dad4f2 100644 --- a/Makefile +++ b/Makefile @@ -28,7 +28,7 @@ shortlog: @hg log -M -r$(RELEASE): --template '- {desc|firstline} (:rev:`{node|short}`)\n' black: - black -l 82 fluidsim + black -l 82 fluidsim scripts bench doc tests: fluidsim-test -v diff --git a/bench/dedalus/ns2d.py b/bench/dedalus/ns2d.py index b997db19..79d78321 100644 --- a/bench/dedalus/ns2d.py +++ b/bench/dedalus/ns2d.py @@ -14,29 +14,31 @@ from dedalus import public as de import time -lx, ly = (1., 1.) +lx, ly = (1.0, 1.0) nx, ny = (512, 512) # Create bases and domain -x_basis = de.Fourier('x', nx, interval=(0, lx), dealias=3/2) -y_basis = de.Fourier('y', ny, interval=(0, ly), dealias=3/2) +x_basis = de.Fourier("x", nx, interval=(0, lx), dealias=3 / 2) +y_basis = de.Fourier("y", ny, interval=(0, ly), dealias=3 / 2) domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64) # Faster algorithm with vorticity? -variables = ['p', 'u', 'v', 'uy', 'vy'] +variables = ["p", "u", "v", "uy", "vy"] problem = de.IVP(domain, variables=variables) Reynolds = 1e4 -problem.parameters['Re'] = Reynolds +problem.parameters["Re"] = Reynolds problem.add_equation( - 'dt(u) + dx(p) - 1/Re*(dx(dx(u)) + dy(uy)) = - u*dx(u) - v*uy') + "dt(u) + dx(p) - 1/Re*(dx(dx(u)) + dy(uy)) = - u*dx(u) - v*uy" +) problem.add_equation( - 'dt(v) + dy(p) - 1/Re*(dx(dx(v)) + dy(vy)) = - u*dx(v) - v*vy') + "dt(v) + dy(p) - 1/Re*(dx(dx(v)) + dy(vy)) = - u*dx(v) - v*vy" +) -problem.add_equation('dx(u) + vy = 0') -problem.add_equation('uy - dy(u) = 0') -problem.add_equation('vy - dy(v) = 0') +problem.add_equation("dx(u) + vy = 0") +problem.add_equation("uy - dy(u) = 0") +problem.add_equation("vy - dy(v) = 0") ts = de.timesteppers.RK443 @@ -44,19 +46,19 @@ x = domain.grid(0) y = domain.grid(1) -u = solver.state['u'] -uy = solver.state['uy'] -v = solver.state['v'] -vy = solver.state['vy'] +u = solver.state["u"] +uy = solver.state["uy"] +v = solver.state["v"] +vy = solver.state["vy"] -u['g'] = np.ones_like(x) -v['g'] = np.ones_like(x) -u.differentiate('y', out=uy) -v.differentiate('y', out=vy) +u["g"] = np.ones_like(x) +v["g"] = np.ones_like(x) +u.differentiate("y", out=uy) +v.differentiate("y", out=vy) dt = 1e-12 -print('Starting loop') +print("Starting loop") start_time = time.time() for it in range(10): @@ -64,4 +66,4 @@ end_time = time.time() -print('Run time: %f' %(end_time-start_time)) +print("Run time: %f" % (end_time - start_time)) diff --git a/bench/dedalus/ns2d_rot.py b/bench/dedalus/ns2d_rot.py index 42d7b166..2050bc06 100644 --- a/bench/dedalus/ns2d_rot.py +++ b/bench/dedalus/ns2d_rot.py @@ -20,32 +20,31 @@ from dedalus import public as de -lx, ly = (1., 1.) +lx, ly = (1.0, 1.0) n = 1024 nx, ny = (n, n) # Create bases and domain -x_basis = de.Fourier('x', nx, interval=(0, lx), dealias=3/2) -y_basis = de.Fourier('y', ny, interval=(0, ly), dealias=3/2) +x_basis = de.Fourier("x", nx, interval=(0, lx), dealias=3 / 2) +y_basis = de.Fourier("y", ny, interval=(0, ly), dealias=3 / 2) domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64) # Stream function-vorticity formulation -variables = ['rot', 'psi', 'u', 'v', 'rotx', 'roty'] +variables = ["rot", "psi", "u", "v", "rotx", "roty"] problem = de.IVP(domain, variables=variables) Reynolds = 1e4 -problem.parameters['Re'] = Reynolds +problem.parameters["Re"] = Reynolds -problem.add_equation( - 'dt(rot) - (1/Re)*(dx(rotx) + dy(roty)) = - u*rotx - v*roty') -problem.add_equation('dx(dx(psi)) + dy(dy(psi)) + rot = 0') +problem.add_equation("dt(rot) - (1/Re)*(dx(rotx) + dy(roty)) = - u*rotx - v*roty") +problem.add_equation("dx(dx(psi)) + dy(dy(psi)) + rot = 0") # with first-order reduction equations... # problem.add_equation('rot + dy(u) - dx(v) = 0') -problem.add_equation('rotx - dx(rot) = 0') -problem.add_equation('roty - dy(rot) = 0') -problem.add_equation('v + dx(psi) = 0') -problem.add_equation('u - dy(psi) = 0') +problem.add_equation("rotx - dx(rot) = 0") +problem.add_equation("roty - dy(rot) = 0") +problem.add_equation("v + dx(psi) = 0") +problem.add_equation("u - dy(psi) = 0") ts = de.timesteppers.RK443 @@ -53,12 +52,12 @@ x = domain.grid(0) y = domain.grid(1) -rot = solver.state['rot'] -psi = solver.state['psi'] -u = solver.state['u'] -v = solver.state['v'] -rotx = solver.state['rotx'] -roty = solver.state['roty'] +rot = solver.state["rot"] +psi = solver.state["psi"] +u = solver.state["u"] +v = solver.state["v"] +rotx = solver.state["rotx"] +roty = solver.state["roty"] # Initial conditions @@ -68,9 +67,9 @@ # rot.differentiate('x', out=rotx) # rot.differentiate('y', out=roty) -rot['g'] = gausspulse(np.sqrt((x - 0.5)**2 + (y - 0.5)**2), fc=1) -rot.differentiate('x', out=rotx) -rot.differentiate('y', out=roty) +rot["g"] = gausspulse(np.sqrt((x - 0.5) ** 2 + (y - 0.5) ** 2), fc=1) +rot.differentiate("x", out=rotx) +rot.differentiate("y", out=roty) # psi['g'] = 10 * y # psi.differentiate('y', out=u) @@ -82,7 +81,7 @@ solver.stop_wall_time = np.inf solver.stop_iteration = 10 -print('Starting main time loop...') +print("Starting main time loop...") start_time = time.time() for it in range(10): @@ -90,4 +89,4 @@ end_time = time.time() -print('Run time for the loop: %f' %(end_time-start_time)) +print("Run time for the loop: %f" % (end_time - start_time)) diff --git a/bench/dedalus/ns2d_rot_faster.py b/bench/dedalus/ns2d_rot_faster.py index f5bd4c16..bee8baea 100644 --- a/bench/dedalus/ns2d_rot_faster.py +++ b/bench/dedalus/ns2d_rot_faster.py @@ -23,35 +23,36 @@ from dedalus import public as de -lx, ly = (1., 1.) +lx, ly = (1.0, 1.0) -nx = 512*2 +nx = 512 * 2 -coef_dealias = 2/3 +coef_dealias = 2 / 3 n = int(coef_dealias * nx) -dealias = nx/n +dealias = nx / n nx, ny = (n, n) # Create bases and domain -x_basis = de.Fourier('x', nx, interval=(0, lx), dealias=dealias) -y_basis = de.Fourier('y', ny, interval=(0, ly), dealias=dealias) +x_basis = de.Fourier("x", nx, interval=(0, lx), dealias=dealias) +y_basis = de.Fourier("y", ny, interval=(0, ly), dealias=dealias) domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64) # Stream function-vorticity formulation -variables = ['psi'] +variables = ["psi"] problem = de.IVP(domain, variables=variables) Reynolds = 1e4 -problem.parameters['Re'] = Reynolds -problem.substitutions['u'] = "dy(psi)" -problem.substitutions['v'] = "-dx(psi)" -problem.substitutions['rot'] = "- dx(dx(psi)) - dy(dy(psi))" -problem.substitutions['rotx'] = "dx(rot)" -problem.substitutions['roty'] = "dy(rot)" +problem.parameters["Re"] = Reynolds +problem.substitutions["u"] = "dy(psi)" +problem.substitutions["v"] = "-dx(psi)" +problem.substitutions["rot"] = "- dx(dx(psi)) - dy(dy(psi))" +problem.substitutions["rotx"] = "dx(rot)" +problem.substitutions["roty"] = "dy(rot)" problem.add_equation( - 'dt(rot) - (1/Re)*(dx(rotx) + dy(roty)) = - u*rotx - v*roty', - condition="(nx != 0) or (ny != 0)") + "dt(rot) - (1/Re)*(dx(rotx) + dy(roty)) = - u*rotx - v*roty", + condition="(nx != 0) or (ny != 0)", +) problem.add_equation("psi = 0", condition="(nx == 0) and (ny == 0)") # with first-order reduction equations... @@ -65,7 +66,7 @@ x = domain.grid(0) y = domain.grid(1) -psi = solver.state['psi'] +psi = solver.state["psi"] # Initial conditions @@ -76,11 +77,11 @@ # rot.differentiate('y', out=roty) rot = domain.new_field() -rot['g'] = gausspulse(np.sqrt((x - 0.5)**2 + (y - 0.5)**2), fc=1) +rot["g"] = gausspulse(np.sqrt((x - 0.5) ** 2 + (y - 0.5) ** 2), fc=1) kx = domain.elements(0) ky = domain.elements(1) -k2 = kx**2 + ky**2 -psi['c'][k2 != 0] = rot['c'][k2 != 0] / k2[k2 != 0] +k2 = kx ** 2 + ky ** 2 +psi["c"][k2 != 0] = rot["c"][k2 != 0] / k2[k2 != 0] # psi['g'] = 10 * y # psi.differentiate('y', out=u) @@ -92,16 +93,16 @@ solver.stop_wall_time = np.inf solver.stop_iteration = 10 -print('Starting startup loop...') +print("Starting startup loop...") start_time = time.time() for it in range(10): solver.step(dt) end_time = time.time() -print('Run time for startup loop: %f' %(end_time-start_time)) +print("Run time for startup loop: %f" % (end_time - start_time)) -print('Starting main time loop...') +print("Starting main time loop...") start_time = time.time() for it in range(10): solver.step(dt) end_time = time.time() -print('Run time for main loop: %f' %(end_time-start_time)) +print("Run time for main loop: %f" % (end_time - start_time)) diff --git a/bench/dedalus/ns2dstrat.py b/bench/dedalus/ns2dstrat.py index 6e3a9750..fbba7481 100644 --- a/bench/dedalus/ns2dstrat.py +++ b/bench/dedalus/ns2dstrat.py @@ -15,44 +15,47 @@ from dedalus import public as de import time -lx, lz = (1., 1.) +lx, lz = (1.0, 1.0) nz, nx = (512, 512) # Create bases and domain -x_basis = de.Fourier('x', nx, interval=(0, lx), dealias=3/2) -z_basis = de.Fourier('z', nz, interval=(0, lz), dealias=3/2) +x_basis = de.Fourier("x", nx, interval=(0, lx), dealias=3 / 2) +z_basis = de.Fourier("z", nz, interval=(0, lz), dealias=3 / 2) domain = de.Domain([x_basis, z_basis], grid_dtype=np.float64) # Faster algorithm with vorticity? -variables = ['p', 'b', 'u', 'w', 'uz', 'wz', 'bz'] +variables = ["p", "b", "u", "w", "uz", "wz", "bz"] problem = de.IVP(domain, variables=variables) # Non-dimensional parameters Reynolds = 1e4 Froude_horiz = 5e-1 -Schmidt = 1. -Aspect_ratio = 1. +Schmidt = 1.0 +Aspect_ratio = 1.0 -problem.parameters['Re'] = Reynolds -problem.parameters['Fh'] = Froude_horiz -problem.parameters['Sc'] = Schmidt -problem.parameters['alpha'] = Aspect_ratio +problem.parameters["Re"] = Reynolds +problem.parameters["Fh"] = Froude_horiz +problem.parameters["Sc"] = Schmidt +problem.parameters["alpha"] = Aspect_ratio # Non-dimensional NS2D of stratified fluid. problem.add_equation( - 'dt(u) + dx(p) - (1/(Re * alpha**2)) * ((alpha**2) * dx(dx(u)) + dz(uz)) = - u*dx(u) - w*uz') + "dt(u) + dx(p) - (1/(Re * alpha**2)) * ((alpha**2) * dx(dx(u)) + dz(uz)) = - u*dx(u) - w*uz" +) problem.add_equation( - '(Fh**2) * dt(w) + dz(p) + b - (1/(Re * alpha**2)) * ((alpha**2) * dx(dx(w)) + dz(wz)) = (Fh**2) * (- u*dx(w) - w*wz)') + "(Fh**2) * dt(w) + dz(p) + b - (1/(Re * alpha**2)) * ((alpha**2) * dx(dx(w)) + dz(wz)) = (Fh**2) * (- u*dx(w) - w*wz)" +) problem.add_equation( - 'dt(b) - uz - (1/(Re * Sc * alpha**2)) * ((alpha**2) * dx(dx(b)) + dz(bz)) = - u*dx(b) - w*bz') + "dt(b) - uz - (1/(Re * Sc * alpha**2)) * ((alpha**2) * dx(dx(b)) + dz(bz)) = - u*dx(b) - w*bz" +) -problem.add_equation('dx(u) + (Fh**2 / alpha**2) * wz = 0') +problem.add_equation("dx(u) + (Fh**2 / alpha**2) * wz = 0") # with first-order reduction equations... -problem.add_equation('uz - dz(u) = 0') -problem.add_equation('wz - dz(w) = 0') -problem.add_equation('bz - dz(b) = 0') +problem.add_equation("uz - dz(u) = 0") +problem.add_equation("wz - dz(w) = 0") +problem.add_equation("bz - dz(b) = 0") ts = de.timesteppers.RK443 @@ -60,23 +63,23 @@ x = domain.grid(0) z = domain.grid(1) -u = solver.state['u'] -uz = solver.state['uz'] -w = solver.state['w'] -wz = solver.state['wz'] -b = solver.state['b'] -bz = solver.state['bz'] - -u['g'] = np.ones_like(x) -w['g'] = np.ones_like(x) -b['g'] = np.ones_like(x) -u.differentiate('z', out=uz) -w.differentiate('z', out=wz) -b.differentiate('z', out=bz) +u = solver.state["u"] +uz = solver.state["uz"] +w = solver.state["w"] +wz = solver.state["wz"] +b = solver.state["b"] +bz = solver.state["bz"] + +u["g"] = np.ones_like(x) +w["g"] = np.ones_like(x) +b["g"] = np.ones_like(x) +u.differentiate("z", out=uz) +w.differentiate("z", out=wz) +b.differentiate("z", out=bz) dt = 1e-12 -print('Starting loop') +print("Starting loop") start_time = time.time() for it in range(10): @@ -84,4 +87,4 @@ end_time = time.time() -print('Run time: %f' %(end_time-start_time)) +print("Run time: %f" % (end_time - start_time)) diff --git a/bench/dedalus/ns2dstrat_rot.py b/bench/dedalus/ns2dstrat_rot.py index d1bf7866..dbc6ee8a 100644 --- a/bench/dedalus/ns2dstrat_rot.py +++ b/bench/dedalus/ns2dstrat_rot.py @@ -15,44 +15,46 @@ from dedalus import public as de import time -lx, lz = (1., 1.) +lx, lz = (1.0, 1.0) nz, nx = (512, 512) # Create bases and domain -x_basis = de.Fourier('x', nx, interval=(0, lx), dealias=3/2) -z_basis = de.Fourier('z', nz, interval=(0, lz), dealias=3/2) +x_basis = de.Fourier("x", nx, interval=(0, lx), dealias=3 / 2) +z_basis = de.Fourier("z", nz, interval=(0, lz), dealias=3 / 2) domain = de.Domain([x_basis, z_basis], grid_dtype=np.float64) # Faster algorithm with vorticity? -variables = ['xi', 'b', 'u', 'w', 'uz', 'wz', 'xiz', 'bz'] +variables = ["xi", "b", "u", "w", "uz", "wz", "xiz", "bz"] problem = de.IVP(domain, variables=variables) # Non-dimensional parameters Reynolds = 1e4 Froude = 5e-1 -Ratio_omegas = 1. -Schmidt = 1. +Ratio_omegas = 1.0 +Schmidt = 1.0 -problem.parameters['Re'] = Reynolds -problem.parameters['F'] = Froude -problem.parameters['R'] = Ratio_omegas -problem.parameters['Sc'] = Schmidt +problem.parameters["Re"] = Reynolds +problem.parameters["F"] = Froude +problem.parameters["R"] = Ratio_omegas +problem.parameters["Sc"] = Schmidt # Non-dimensional NS2D of stratified fluid. problem.add_equation( - 'dt(xi) + (R**2 / (F * (1 - F**2)**(1/2.))) * dx(b) - (F**2 / Re) * dx(dx(xi)) - ((1 - F**2) / Re) * dz(xiz) = - u * dx(xi) - w * xiz') + "dt(xi) + (R**2 / (F * (1 - F**2)**(1/2.))) * dx(b) - (F**2 / Re) * dx(dx(xi)) - ((1 - F**2) / Re) * dz(xiz) = - u * dx(xi) - w * xiz" +) problem.add_equation( - 'dt(b) + uz - (F**2/(Re * Sc)) * dx(dx(b) - ((1 - F**2)/(Re * Sc))*dz(bz)) = - u*dx(b) - w*bz') + "dt(b) + uz - (F**2/(Re * Sc)) * dx(dx(b) - ((1 - F**2)/(Re * Sc))*dz(bz)) = - u*dx(b) - w*bz" +) -problem.add_equation('dx(u) + wz = 0') -problem.add_equation('xi - dz(u) + dx(w) = 0') +problem.add_equation("dx(u) + wz = 0") +problem.add_equation("xi - dz(u) + dx(w) = 0") # with first-order reduction equations... -problem.add_equation('uz - dz(u) = 0') -problem.add_equation('wz - dz(w) = 0') -problem.add_equation('bz - dz(b) = 0') -problem.add_equation('xiz - dz(xi) = 0') +problem.add_equation("uz - dz(u) = 0") +problem.add_equation("wz - dz(w) = 0") +problem.add_equation("bz - dz(b) = 0") +problem.add_equation("xiz - dz(xi) = 0") ts = de.timesteppers.RK443 @@ -60,27 +62,27 @@ x = domain.grid(0) z = domain.grid(1) -xi = solver.state['xi'] -xiz = solver.state['xiz'] -u = solver.state['u'] -uz = solver.state['uz'] -w = solver.state['w'] -wz = solver.state['wz'] -b = solver.state['b'] -bz = solver.state['bz'] - -u['g'] = np.ones_like(x) -w['g'] = np.ones_like(x) -b['g'] = np.ones_like(x) -xi['g'] = np.ones_like(x) -u.differentiate('z', out=uz) -w.differentiate('z', out=wz) -b.differentiate('z', out=bz) -xi.differentiate('z', out=xiz) +xi = solver.state["xi"] +xiz = solver.state["xiz"] +u = solver.state["u"] +uz = solver.state["uz"] +w = solver.state["w"] +wz = solver.state["wz"] +b = solver.state["b"] +bz = solver.state["bz"] + +u["g"] = np.ones_like(x) +w["g"] = np.ones_like(x) +b["g"] = np.ones_like(x) +xi["g"] = np.ones_like(x) +u.differentiate("z", out=uz) +w.differentiate("z", out=wz) +b.differentiate("z", out=bz) +xi.differentiate("z", out=xiz) dt = 1e-12 -print('Starting loop') +print("Starting loop") start_time = time.time() for it in range(10): @@ -88,4 +90,4 @@ end_time = time.time() -print('Run time: %f' %(end_time-start_time)) +print("Run time: %f" % (end_time - start_time)) diff --git a/bench/old/compare_fluidfft.py b/bench/old/compare_fluidfft.py index 7e8d1628..db7492cf 100644 --- a/bench/old/compare_fluidfft.py +++ b/bench/old/compare_fluidfft.py @@ -15,33 +15,35 @@ def modif_params(params, old=False): - params.short_name_type_run = 'profile' + params.short_name_type_run = "profile" - nh = 512*2 + nh = 512 * 2 params.oper.nx = nh params.oper.ny = nh - lh = 6. + lh = 6.0 params.oper.Lx = lh params.oper.Ly = lh - params.oper.coef_dealiasing = 2./3 + params.oper.coef_dealiasing = 2.0 / 3 params.FORCING = True - params.forcing.type = 'tcrandom' + params.forcing.type = "tcrandom" params.forcing.nkmax_forcing = 5 params.forcing.nkmin_forcing = 4 - params.forcing.forcing_rate = 1. + params.forcing.forcing_rate = 1.0 - delta_x = lh/nh - params.nu_8 = 2.*10e-1*params.forcing.forcing_rate**(1./3)*delta_x**8 + delta_x = lh / nh + params.nu_8 = ( + 2.0 * 10e-1 * params.forcing.forcing_rate ** (1.0 / 3) * delta_x ** 8 + ) try: - params.f = 1. - params.c2 = 200. + params.f = 1.0 + params.c2 = 200.0 except (KeyError, AttributeError): pass - params.time_stepping.deltat0 = 1.e-4 + params.time_stepping.deltat0 = 1.0e-4 params.time_stepping.USE_CFL = False params.time_stepping.it_end = 10 @@ -57,20 +59,24 @@ def modif_params(params, old=False): # params.output.periods_save.increments = 0. if not old: - params.short_name_type_run = 'profile2' + params.short_name_type_run = "profile2" # params.oper.type_fft = 'fft2d.mpi_with_fftw1d' else: - params.oper.type_fft = 'FFTWCY' + params.oper.type_fft = "FFTWCY" params = SimulOperCython.create_default_params() -modif_params(params, 'old') +modif_params(params, "old") sim_oper_cython = SimulOperCython(params) t_start0 = time() -cProfile.runctx('sim_oper_cython.time_stepping.start()', - globals(), locals(), 'profile_oper_cython.pstats') +cProfile.runctx( + "sim_oper_cython.time_stepping.start()", + globals(), + locals(), + "profile_oper_cython.pstats", +) t_end0 = time() @@ -79,23 +85,28 @@ def modif_params(params, old=False): sim = Simul(params) t_start1 = time() -cProfile.runctx('sim.time_stepping.start()', - globals(), locals(), 'profile.pstats') +cProfile.runctx( + "sim.time_stepping.start()", globals(), locals(), "profile.pstats" +) t_end1 = time() if sim.oper.rank == 0: - s = pstats.Stats('profile_oper_cython.pstats') - s.strip_dirs().sort_stats('time').print_stats(10) + s = pstats.Stats("profile_oper_cython.pstats") + s.strip_dirs().sort_stats("time").print_stats(10) - s = pstats.Stats('profile.pstats') - s.strip_dirs().sort_stats('time').print_stats(10) + s = pstats.Stats("profile.pstats") + s.strip_dirs().sort_stats("time").print_stats(10) - print('elapsed times: {:.3f} and {:.3f}'.format( - t_end0 - t_start0, t_end1 - t_start1)) + print( + "elapsed times: {:.3f} and {:.3f}".format( + t_end0 - t_start0, t_end1 - t_start1 + ) + ) print( - 'with gprof2dot and graphviz (command dot):\n' - 'gprof2dot -f pstats profile.pstats | dot -Tpng -o profile.png\n' - 'gprof2dot -f pstats profile_oper_cython.pstats | ' - 'dot -Tpng -o profile_oper_cython.png') + "with gprof2dot and graphviz (command dot):\n" + "gprof2dot -f pstats profile.pstats | dot -Tpng -o profile.png\n" + "gprof2dot -f pstats profile_oper_cython.pstats | " + "dot -Tpng -o profile_oper_cython.png" + ) diff --git a/bench/old/plot_bench.py b/bench/old/plot_bench.py index 4c434ece..78965c81 100644 --- a/bench/old/plot_bench.py +++ b/bench/old/plot_bench.py @@ -1,38 +1,36 @@ - from glob import glob import json import pandas as pd import matplotlib.pyplot as plt -key_solver = 'ns2d' -hostname = 'cl7' +key_solver = "ns2d" +hostname = "cl7" dicts = [] -for path in glob('results_bench/*'): +for path in glob("results_bench/*"): with open(path) as file: d = json.load(file) - if not d['hostname'].startswith(hostname): + if not d["hostname"].startswith(hostname): continue dicts.append(d) - df = pd.DataFrame(dicts) def plot_bench(nb_proc0=1): - df0 = df.loc[df['nb_proc'] == nb_proc0] - t_elapsed0 = df0['t_elapsed'].mean() - times = df['t_elapsed'] + df0 = df.loc[df["nb_proc"] == nb_proc0] + t_elapsed0 = df0["t_elapsed"].mean() + times = df["t_elapsed"] plt.figure() ax = plt.subplot() - ax.plot(df['nb_proc'], nb_proc0*t_elapsed0/times, 'xr') - tmp = [nb_proc0, df['nb_proc'].max()] - ax.plot(tmp, tmp, 'b-') - ax.set_title('speed up') + ax.plot(df["nb_proc"], nb_proc0 * t_elapsed0 / times, "xr") + tmp = [nb_proc0, df["nb_proc"].max()] + ax.plot(tmp, tmp, "b-") + ax.set_title("speed up") plot_bench(nb_proc0=1) diff --git a/bench/old/simul_bench.py b/bench/old/simul_bench.py index 8c393a2c..e3d25b29 100644 --- a/bench/old/simul_bench.py +++ b/bench/old/simul_bench.py @@ -7,17 +7,17 @@ from fluidsim import import_module_solver_from_key from util_bench import modif_params_profile2d, modif_params_profile3d, bench -key = 'ns2d' +key = "ns2d" solver = import_module_solver_from_key(key) params = solver.Simul.create_default_params() -if '3d' in key: +if "3d" in key: modif_params_profile3d(params) else: modif_params_profile2d(params) sim = solver.Simul(params) -if __name__ == '__main__': +if __name__ == "__main__": bench(sim) diff --git a/bench/old/submit_benchs.py b/bench/old/submit_benchs.py index e0a49729..552baa28 100644 --- a/bench/old/submit_benchs.py +++ b/bench/old/submit_benchs.py @@ -1,5 +1,5 @@ - from fluiddyn.clusters.legi import Calcul8 as Cluster + cluster = Cluster() # cluster.submit_script( @@ -9,7 +9,11 @@ for nb_mpi_processes in [2, 4, 8, 16]: cluster.submit_script( - 'simul_bench.py', name_run='fluidsim_bench', walltime='1:00:00', + "simul_bench.py", + name_run="fluidsim_bench", + walltime="1:00:00", omp_num_threads=1, nb_cores_per_node=nb_mpi_processes, - nb_mpi_processes=nb_mpi_processes, ask=False) + nb_mpi_processes=nb_mpi_processes, + ask=False, + ) diff --git a/bench/profile_time_stepping/init.py b/bench/profile_time_stepping/init.py index 58610e86..3f2c537e 100644 --- a/bench/profile_time_stepping/init.py +++ b/bench/profile_time_stepping/init.py @@ -1,14 +1,13 @@ - from fluidsim.solvers.ns2d.solver import Simul from fluidsim.solvers.sw1l.solver import Simul params = Simul.create_default_params() -params.output.sub_directory = 'bench_time_stepping' +params.output.sub_directory = "bench_time_stepping" params.oper.nx = params.oper.ny = nh = 1024 -params.nu_8 = 1. +params.nu_8 = 1.0 params.time_stepping.it_end = 10 params.time_stepping.USE_T_END = False @@ -16,10 +15,10 @@ params.time_stepping.deltat0 = 1e-16 params.time_stepping.type_time_scheme = "RK2" -params.init_fields.type = 'dipole' +params.init_fields.type = "dipole" params.forcing.enable = True -params.forcing.type = 'proportional' +params.forcing.type = "proportional" params.output.periods_print.print_stdout = 0.25 diff --git a/bench/profiling/simul_profile2d.py b/bench/profiling/simul_profile2d.py index e5ae2259..37d1c61d 100644 --- a/bench/profiling/simul_profile2d.py +++ b/bench/profiling/simul_profile2d.py @@ -9,7 +9,7 @@ from fluidsim import import_module_solver_from_key from util_bench import profile, modif_params_profile2d -key = 'ns2d' +key = "ns2d" solver = import_module_solver_from_key(key) @@ -18,5 +18,5 @@ sim = solver.Simul(params) -if __name__ == '__main__': +if __name__ == "__main__": profile(sim) diff --git a/bench/profiling/simul_profile_ns2d.py b/bench/profiling/simul_profile_ns2d.py index 5843e3a6..e58ff357 100644 --- a/bench/profiling/simul_profile_ns2d.py +++ b/bench/profiling/simul_profile_ns2d.py @@ -16,5 +16,5 @@ sim = solver.Simul(params) -if __name__ == '__main__': +if __name__ == "__main__": profile(sim) diff --git a/bench/profiling/simul_profile_ns2dstrat.py b/bench/profiling/simul_profile_ns2dstrat.py index 08234046..207694eb 100644 --- a/bench/profiling/simul_profile_ns2dstrat.py +++ b/bench/profiling/simul_profile_ns2dstrat.py @@ -16,5 +16,5 @@ sim = solver.Simul(params) -if __name__ == '__main__': +if __name__ == "__main__": profile(sim) diff --git a/bench/profiling/simul_profile_ns3d.py b/bench/profiling/simul_profile_ns3d.py index c7d3ea36..a9ba096c 100644 --- a/bench/profiling/simul_profile_ns3d.py +++ b/bench/profiling/simul_profile_ns3d.py @@ -13,11 +13,10 @@ sim = solver.Simul(params) -if __name__ == '__main__': +if __name__ == "__main__": profile(sim, nb_dim=3) # params.oper.type_fft = 'fluidfft.fft3d.with_fftw3d' # params.oper.type_fft = 'fluidfft.fft3d.with_cufft' # params.oper.type_fft = 'fluidfft.fft3d.mpi_with_fftwmpi3d' - diff --git a/bench/profiling/simul_profile_plate2d.py b/bench/profiling/simul_profile_plate2d.py index fce23969..3fb1049f 100644 --- a/bench/profiling/simul_profile_plate2d.py +++ b/bench/profiling/simul_profile_plate2d.py @@ -17,5 +17,5 @@ sim = solver.Simul(params) -if __name__ == '__main__': +if __name__ == "__main__": profile(sim) diff --git a/bench/profiling/simul_profile_sw1l.py b/bench/profiling/simul_profile_sw1l.py index 5976cac9..508cc505 100644 --- a/bench/profiling/simul_profile_sw1l.py +++ b/bench/profiling/simul_profile_sw1l.py @@ -16,5 +16,5 @@ sim = solver.Simul(params) -if __name__ == '__main__': +if __name__ == "__main__": profile(sim) diff --git a/bench/profiling/util_bench.py b/bench/profiling/util_bench.py index 6a586bde..f96d4715 100644 --- a/bench/profiling/util_bench.py +++ b/bench/profiling/util_bench.py @@ -1,4 +1,3 @@ - import os from time import time import pstats @@ -10,8 +9,8 @@ from fluiddyn.util import mpi -def modif_params_profile2d(params, nh=3**2*2**7): - params.short_name_type_run = 'profile' +def modif_params_profile2d(params, nh=3 ** 2 * 2 ** 7): + params.short_name_type_run = "profile" params.oper.nx = nh params.oper.ny = nh @@ -20,24 +19,24 @@ def modif_params_profile2d(params, nh=3**2*2**7): # params.oper.type_fft = 'fft2d.with_cufft' params.forcing.enable = True - params.forcing.type = 'tcrandom' + params.forcing.type = "tcrandom" params.forcing.nkmax_forcing = 5 params.forcing.nkmin_forcing = 4 - params.forcing.forcing_rate = 1. + params.forcing.forcing_rate = 1.0 - params.nu_8 = 1. + params.nu_8 = 1.0 try: - params.f = 1. - params.c2 = 200. + params.f = 1.0 + params.c2 = 200.0 except AttributeError: pass try: - params.N = 1. + params.N = 1.0 except AttributeError: pass - params.time_stepping.deltat0 = 1.e-4 + params.time_stepping.deltat0 = 1.0e-4 params.time_stepping.USE_CFL = False params.time_stepping.it_end = 10 params.time_stepping.USE_T_END = False @@ -47,14 +46,14 @@ def modif_params_profile2d(params, nh=3**2*2**7): def modif_params_profile3d(params, nh=256, nz=32): - params.short_name_type_run = 'profile' + params.short_name_type_run = "profile" params.oper.nx = nh params.oper.ny = nh params.oper.nz = nz # params.oper.type_fft = 'fft3d.with_fftw3d' - params.oper.type_fft = 'fft3d.with_cufft' + params.oper.type_fft = "fft3d.with_cufft" # params.forcing.enable = False # params.forcing.type = 'tcrandom' @@ -62,19 +61,19 @@ def modif_params_profile3d(params, nh=256, nz=32): # params.forcing.nkmin_forcing = 4 # params.forcing.forcing_rate = 1. - params.nu_8 = 1. + params.nu_8 = 1.0 try: - params.f = 1. - params.c2 = 200. + params.f = 1.0 + params.c2 = 200.0 except AttributeError: pass try: - params.N = 1. + params.N = 1.0 except AttributeError: pass - params.time_stepping.deltat0 = 1.e-4 + params.time_stepping.deltat0 = 1.0e-4 params.time_stepping.USE_CFL = False params.time_stepping.it_end = 10 params.time_stepping.USE_T_END = False @@ -87,13 +86,14 @@ def profile(sim, nb_dim=2): t0 = time() - cProfile.runctx('sim.time_stepping.start()', - globals(), locals(), 'profile.pstats') + cProfile.runctx( + "sim.time_stepping.start()", globals(), locals(), "profile.pstats" + ) t_end = time() if sim.oper.rank == 0: - s = pstats.Stats('profile.pstats') + s = pstats.Stats("profile.pstats") # s.strip_dirs().sort_stats('time').print_stats(16) - s.sort_stats('time').print_stats(12) + s.sort_stats("time").print_stats(12) if nb_dim == 2: times = print_analysis(s) @@ -102,11 +102,12 @@ def profile(sim, nb_dim=2): else: raise NotImplementedError - print('\nelapsed time = {:.3f} s'.format(t_end - t0)) + print("\nelapsed time = {:.3f} s".format(t_end - t0)) print( - '\nwith gprof2dot and graphviz (command dot):\n' - 'gprof2dot -f pstats profile.pstats | dot -Tpng -o profile.png') + "\nwith gprof2dot and graphviz (command dot):\n" + "gprof2dot -f pstats profile.pstats | dot -Tpng -o profile.png" + ) def bench(sim): @@ -119,80 +120,93 @@ def bench(sim): if sim.oper.rank != 0: return - path_results = 'results_bench' + path_results = "results_bench" if not os.path.exists(path_results): os.mkdir(path_results) key_solver = sim.info_solver.short_name.lower() pid = os.getpid() nfile = ( - 'result_bench_' + key_solver + '_' + - sim.oper.produce_str_describing_grid() + - '_' + t_as_str + f'_{pid}' + '.json') + "result_bench_" + + key_solver + + "_" + + sim.oper.produce_str_describing_grid() + + "_" + + t_as_str + + f"_{pid}" + + ".json" + ) path = os.path.join(path_results, nfile) results = { - 't_elapsed': t_elapsed, - 'key_solver': sim.info_solver.short_name.lower(), - 'nx': sim.oper.nx, - 'ny': sim.oper.ny, - 'nb_proc': mpi.nb_proc, - 'pid': pid, - 'time_as_str': t_as_str, - 'hostname': socket.gethostname()} + "t_elapsed": t_elapsed, + "key_solver": sim.info_solver.short_name.lower(), + "nx": sim.oper.nx, + "ny": sim.oper.ny, + "nb_proc": mpi.nb_proc, + "pid": pid, + "time_as_str": t_as_str, + "hostname": socket.gethostname(), + } try: - results['nz'] = sim.oper.nz + results["nz"] = sim.oper.nz except AttributeError: pass - with open(path, 'w') as file: + with open(path, "w") as file: json.dump(results, file, sort_keys=True) - file.write('\n') + file.write("\n") def print_analysis(s): - total_time = 0. - times = {'fft2d': 0., 'fft_as': 0., 'pythran': 0., '.pyx': 0.} + total_time = 0.0 + times = {"fft2d": 0.0, "fft_as": 0.0, "pythran": 0.0, ".pyx": 0.0} for key, value in s.stats.items(): name = key[2] time = value[2] total_time += time - if name == 'one_time_step_computation': - print('warning: special case one_time_step_computation ' - 'included in .pyx (see explanation in the code)') - times['.pyx'] += time + if name == "one_time_step_computation": + print( + "warning: special case one_time_step_computation " + "included in .pyx (see explanation in the code)" + ) + times[".pyx"] += time for k in times.keys(): if k in name or k in key[0]: - if k == '.pyx': - if 'fft/Sources' in key[0]: + if k == ".pyx": + if "fft/Sources" in key[0]: continue - if 'fft_as_arg' in key[2]: + if "fft_as_arg" in key[2]: continue - if k == 'fft2d': + if k == "fft2d": - if 'util_pythran' in key[2] or \ - 'operators.py' in key[0] or \ - 'fft_as_arg' in key[2]: + if ( + "util_pythran" in key[2] + or "operators.py" in key[0] + or "fft_as_arg" in key[2] + ): continue callers = value[4] time = 0 for kcaller, vcaller in callers.items(): - if 'fft_as_arg' not in kcaller[2] and\ - 'fft_as_arg' not in kcaller[0]: + if ( + "fft_as_arg" not in kcaller[2] + and "fft_as_arg" not in kcaller[0] + ): time += vcaller[2] # print(k, key) # print(value[:100]) # print(time, '\n') - if k == 'fft_as': - if '.pyx' in key[0]: + if k == "fft_as": + if ".pyx" in key[0]: continue # time = value[3] @@ -202,43 +216,52 @@ def print_analysis(s): times[k] += time - print('Analysis (percentage of total time):') + print("Analysis (percentage of total time):") keys = list(times.keys()) keys.sort(key=lambda key: times[key], reverse=True) for k in keys: t = times[k] - print('time {:10s}: {:5.01f} % ({:4.02f} s)'.format( - k, t/total_time*100, t)) - - print('-' * 24 + '\n{:15s} {:5.01f} %'.format( - '', sum([t for t in times.values()])/total_time*100)) + print( + "time {:10s}: {:5.01f} % ({:4.02f} s)".format( + k, t / total_time * 100, t + ) + ) + + print( + "-" * 24 + + "\n{:15s} {:5.01f} %".format( + "", sum([t for t in times.values()]) / total_time * 100 + ) + ) return times def print_analysis3d(s): - total_time = 0. - times = {'fft3d': 0., 'fft_as': 0., 'pythran': 0., '.pyx': 0.} + total_time = 0.0 + times = {"fft3d": 0.0, "fft_as": 0.0, "pythran": 0.0, ".pyx": 0.0} for key, value in s.stats.items(): name = key[2] time = value[2] total_time += time - if name == 'one_time_step_computation': - print('warning: special case one_time_step_computation ' - 'included in .pyx (see explanation in the code)') - times['.pyx'] += time + if name == "one_time_step_computation": + print( + "warning: special case one_time_step_computation " + "included in .pyx (see explanation in the code)" + ) + times[".pyx"] += time for k in times.keys(): if k in name or k in key[0]: - if k == 'fft3d': - if 'pythran' in key[0] or 'pythran' in key[2]: + if k == "fft3d": + if "pythran" in key[0] or "pythran" in key[2]: continue - if 'operators.py' in key[0]: + if "operators.py" in key[0]: continue - if 'as_arg' in key[2]: + if "as_arg" in key[2]: continue # print(k, key) @@ -247,17 +270,24 @@ def print_analysis3d(s): times[k] += time - print('Analysis (percentage of total time):') + print("Analysis (percentage of total time):") keys = list(times.keys()) keys.sort(key=lambda key: times[key], reverse=True) for k in keys: t = times[k] - print('time {:10s}: {:5.01f} % ({:4.02f} s)'.format( - k, t/total_time*100, t)) - - print('-' * 24 + '\n{:15s} {:5.01f} %'.format( - '', sum([t for t in times.values()])/total_time*100)) + print( + "time {:10s}: {:5.01f} % ({:4.02f} s)".format( + k, t / total_time * 100, t + ) + ) + + print( + "-" * 24 + + "\n{:15s} {:5.01f} %".format( + "", sum([t for t in times.values()]) / total_time * 100 + ) + ) return times diff --git a/bench/spectraldns/bench_2d.py b/bench/spectraldns/bench_2d.py index 1ebea1e0..d2717ec8 100644 --- a/bench/spectraldns/bench_2d.py +++ b/bench/spectraldns/bench_2d.py @@ -57,25 +57,25 @@ def update(context): plt.pause(1e-6) -if __name__ == '__main__': +if __name__ == "__main__": Re = 1e4 U = 2 ** 0.5 - L = 1. + L = 1.0 config.update( { - 'nu': U * L / Re, - 'dt': 1e-12, - 'T': 1e-11, # Should run 10 iterations - 'write_result': 100, - 'L': [L, L], - 'M': [10, 10] # Mesh size is pow(2, M[i]) in direction i + "nu": U * L / Re, + "dt": 1e-12, + "T": 1e-11, # Should run 10 iterations + "write_result": 100, + "L": [L, L], + "M": [10, 10] # Mesh size is pow(2, M[i]) in direction i # 2**9 == 512 - }, 'doublyperiodic' + }, + "doublyperiodic", ) # required to allow overloading through commandline - config.doublyperiodic.add_argument( - "--plot_result", type=int, default=10) + config.doublyperiodic.add_argument("--plot_result", type=int, default=10) if plt is None: sol = get_solver(mesh="doublyperiodic") else: @@ -88,4 +88,4 @@ def update(context): start_time = time.time() solve(sol, context) end_time = time.time() - print('Run time: %f' % (end_time - start_time)) + print("Run time: %f" % (end_time - start_time)) diff --git a/bench/spectraldns/bench_3d.py b/bench/spectraldns/bench_3d.py index 2098f290..d7a692c7 100644 --- a/bench/spectraldns/bench_3d.py +++ b/bench/spectraldns/bench_3d.py @@ -19,6 +19,7 @@ import time import numpy as np from numpy import zeros, sum, float64, sin, cos, prod, asscalar + # from numpy.linalg import norm from spectralDNS import config, get_solver, solve @@ -27,7 +28,7 @@ def initialize(solver, context): - if 'NS' in config.params.solver: + if "NS" in config.params.solver: initialize1(solver, context) else: initialize2(solver, context) @@ -58,8 +59,11 @@ def initialize2(solver, context): def energy_fourier(comm, a): N = config.params.N - result = 2 * sum(abs(a[..., 1:-1])**2) + \ - sum(abs(a[..., 0])**2) + sum(abs(a[..., -1])**2) + result = ( + 2 * sum(abs(a[..., 1:-1]) ** 2) + + sum(abs(a[..., 0]) ** 2) + + sum(abs(a[..., -1]) ** 2) + ) result = comm.allreduce(result) return result @@ -76,40 +80,51 @@ def update(context): params = config.params solver = config.solver - if (params.tstep % params.compute_energy == 0 or - params.tstep % params.plot_step == 0 and params.plot_step > 0): + if ( + params.tstep % params.compute_energy == 0 + or params.tstep % params.plot_step == 0 + and params.plot_step > 0 + ): U = solver.get_velocity(**c) curl = solver.get_curl(**c) - if params.solver == 'NS': + if params.solver == "NS": P = solver.get_pressure(**c) if plt is not None: - if params.tstep % params.plot_step == 0 and solver.rank == 0 and params.plot_step > 0: + if ( + params.tstep % params.plot_step == 0 + and solver.rank == 0 + and params.plot_step > 0 + ): if im1 is None: plt.figure() - im1 = plt.contourf(c.X[1][:, :, 0], c.X[0] - [:, :, 0], U[0, :, :, 10], 100) + im1 = plt.contourf( + c.X[1][:, :, 0], c.X[0][:, :, 0], U[0, :, :, 10], 100 + ) plt.colorbar(im1) plt.draw() globals().update(im1=im1) else: im1.ax.clear() - im1.ax.contourf(c.X[1][:, :, 0], c.X[0] - [:, :, 0], U[0, :, :, 10], 100) + im1.ax.contourf( + c.X[1][:, :, 0], c.X[0][:, :, 0], U[0, :, :, 10], 100 + ) im1.autoscale() plt.pause(1e-6) if params.tstep % params.compute_energy == 0: dx, L = params.dx, params.L ww = solver.comm.reduce( - sum(curl.astype(float64) * curl.astype(float64)) / prod(params.N) / 2) + sum(curl.astype(float64) * curl.astype(float64)) / prod(params.N) / 2 + ) # Compute energy with double precision kk = solver.comm.reduce( - sum(U.astype(float64) * U.astype(float64)) / prod(params.N) / 2) - if 'shenfun' in params.solver: + sum(U.astype(float64) * U.astype(float64)) / prod(params.N) / 2 + ) + if "shenfun" in params.solver: ww2 = energy_fourier(solver.comm, c.U_hat) / 2 else: - ww2 = energy_fourier(solver.comm, c.U_hat) / prod(params.N)**2 / 2 + ww2 = energy_fourier(solver.comm, c.U_hat) / prod(params.N) ** 2 / 2 kold[0] = kk if solver.rank == 0: @@ -125,11 +140,13 @@ def regression_test(context): U = solver.get_velocity(**context) curl = solver.get_curl(**context) w = solver.comm.reduce( - sum(curl.astype(float64) * curl.astype(float64)) / prod(params.N) / 2) + sum(curl.astype(float64) * curl.astype(float64)) / prod(params.N) / 2 + ) # Compute energy with double precision k = solver.comm.reduce( - sum(U.astype(float64) * U.astype(float64)) / prod(params.N) / 2) - config.solver.MemoryUsage('End', solver.comm) + sum(U.astype(float64) * U.astype(float64)) / prod(params.N) / 2 + ) + config.solver.MemoryUsage("End", solver.comm) if solver.rank == 0: assert round(asscalar(w) - 0.375249930801, params.ntol) == 0 assert round(asscalar(k) - 0.124953117517, params.ntol) == 0 @@ -137,35 +154,37 @@ def regression_test(context): if __name__ == "__main__": Re = 1e4 - U = 2 ** (1. / 3) - L = 1. + U = 2 ** (1.0 / 3) + L = 1.0 dt = 1e-12 config.update( { - 'nu': U * L / Re, # Viscosity - 'dt': dt, # Time step - 'T': 11*dt, # End time - 'L': [L, L, L], - 'M': [7, 7, 7], # Mesh size is pow(2, M[i]) in direction i + "nu": U * L / Re, # Viscosity + "dt": dt, # Time step + "T": 11 * dt, # End time + "L": [L, L, L], + "M": [7, 7, 7], # Mesh size is pow(2, M[i]) in direction i #'planner_effort': {'fft': 'FFTW_EXHAUSTIVE'}, #'decomposition': 'pencil', #'P1': 2 - }, "triplyperiodic" + }, + "triplyperiodic", ) config.triplyperiodic.add_argument("--compute_energy", type=int, default=2) config.triplyperiodic.add_argument("--plot_step", type=int, default=2) if plt is None: sol = get_solver(mesh="triplyperiodic") else: - sol = get_solver(update=update, regression_test=regression_test, - mesh="triplyperiodic") + sol = get_solver( + update=update, regression_test=regression_test, mesh="triplyperiodic" + ) context = sol.get_context() # Add curl to the stored results. For this we need to update the update_components # method used by the HDF5Writer class to compute the real fields that are stored WRITE_HDF5 = False - if config.params.solver == 'NS' and WRITE_HDF5: + if config.params.solver == "NS" and WRITE_HDF5: context.hdf5file.fname = "NS8.h5" context.hdf5file.components["curlx"] = context.curl[0] context.hdf5file.components["curly"] = context.curl[1] @@ -185,7 +204,7 @@ def update_components(**context): start_time = time.time() solve(sol, context) end_time = time.time() - print('Run time: %f' % (end_time - start_time)) + print("Run time: %f" % (end_time - start_time)) - #context.hdf5file._init_h5file(config.params, **context) + # context.hdf5file._init_h5file(config.params, **context) # context.hdf5file.f.close() diff --git a/bench/spectraldns/profile_2d.py b/bench/spectraldns/profile_2d.py index 2fd2ae5d..c0b543ec 100644 --- a/bench/spectraldns/profile_2d.py +++ b/bench/spectraldns/profile_2d.py @@ -59,26 +59,26 @@ def update(context): plt.pause(1e-6) -if __name__ == '__main__': +if __name__ == "__main__": Re = 1e4 U = 2 ** 0.5 - L = 1. + L = 1.0 dt = 1e-12 config.update( { - 'nu': U * L / Re, - 'dt': dt, - 'T': 11 * dt, # Should run 10 iterations - 'write_result': 100, - 'L': [L, L], - 'M': [10, 10] # Mesh size is pow(2, M[i]) in direction i + "nu": U * L / Re, + "dt": dt, + "T": 11 * dt, # Should run 10 iterations + "write_result": 100, + "L": [L, L], + "M": [10, 10] # Mesh size is pow(2, M[i]) in direction i # 2**9 == 512 - }, 'doublyperiodic' + }, + "doublyperiodic", ) # required to allow overloading through commandline - config.doublyperiodic.add_argument( - "--plot_result", type=int, default=10) + config.doublyperiodic.add_argument("--plot_result", type=int, default=10) if plt is None: sol = get_solver(mesh="doublyperiodic") else: @@ -89,13 +89,14 @@ def update(context): # Double check benchmark walltime start_time = time.time() - cProfile.runctx('solve(sol, context)', - globals(), locals(), 'profile.pstats') + cProfile.runctx("solve(sol, context)", globals(), locals(), "profile.pstats") end_time = time.time() - print('Run time: %f' % (end_time - start_time)) + print("Run time: %f" % (end_time - start_time)) - s = pstats.Stats('profile.pstats') - s.sort_stats('time').print_stats(12) + s = pstats.Stats("profile.pstats") + s.sort_stats("time").print_stats(12) - print('you can run:\n' - 'gprof2dot -f pstats profile.pstats | dot -Tpng -o profile.png') + print( + "you can run:\n" + "gprof2dot -f pstats profile.pstats | dot -Tpng -o profile.png" + ) diff --git a/bench/spectraldns/profile_3d.py b/bench/spectraldns/profile_3d.py index 87377401..c89473d1 100644 --- a/bench/spectraldns/profile_3d.py +++ b/bench/spectraldns/profile_3d.py @@ -28,7 +28,7 @@ def initialize(solver, context): - if 'NS' in config.params.solver: + if "NS" in config.params.solver: initialize1(solver, context) else: initialize2(solver, context) @@ -59,9 +59,11 @@ def initialize2(solver, context): def energy_fourier(comm, a): # N = config.params.N - result = 2 * sum(abs(a[..., 1:-1])**2) \ - + sum(abs(a[..., 0])**2) \ - + sum(abs(a[..., -1])**2) + result = ( + 2 * sum(abs(a[..., 1:-1]) ** 2) + + sum(abs(a[..., 0]) ** 2) + + sum(abs(a[..., -1]) ** 2) + ) result = comm.allreduce(result) return result @@ -78,40 +80,51 @@ def update(context): params = config.params solver = config.solver - if (params.tstep % params.compute_energy == 0 or - params.tstep % params.plot_step == 0 and params.plot_step > 0): + if ( + params.tstep % params.compute_energy == 0 + or params.tstep % params.plot_step == 0 + and params.plot_step > 0 + ): U = solver.get_velocity(**c) curl = solver.get_curl(**c) # if params.solver == 'NS': # P = solver.get_pressure(**c) if plt is not None: - if params.tstep % params.plot_step == 0 and solver.rank == 0 and params.plot_step > 0: + if ( + params.tstep % params.plot_step == 0 + and solver.rank == 0 + and params.plot_step > 0 + ): if im1 is None: plt.figure() - im1 = plt.contourf(c.X[1][:, :, 0], c.X[0] - [:, :, 0], U[0, :, :, 10], 100) + im1 = plt.contourf( + c.X[1][:, :, 0], c.X[0][:, :, 0], U[0, :, :, 10], 100 + ) plt.colorbar(im1) plt.draw() globals().update(im1=im1) else: im1.ax.clear() - im1.ax.contourf(c.X[1][:, :, 0], c.X[0] - [:, :, 0], U[0, :, :, 10], 100) + im1.ax.contourf( + c.X[1][:, :, 0], c.X[0][:, :, 0], U[0, :, :, 10], 100 + ) im1.autoscale() plt.pause(1e-6) if params.tstep % params.compute_energy == 0: # dx, L = params.dx, params.L ww = solver.comm.reduce( - sum(curl.astype(float64) * curl.astype(float64)) / prod(params.N) / 2) + sum(curl.astype(float64) * curl.astype(float64)) / prod(params.N) / 2 + ) # Compute energy with double precision kk = solver.comm.reduce( - sum(U.astype(float64) * U.astype(float64)) / prod(params.N) / 2) - if 'shenfun' in params.solver: + sum(U.astype(float64) * U.astype(float64)) / prod(params.N) / 2 + ) + if "shenfun" in params.solver: ww2 = energy_fourier(solver.comm, c.U_hat) / 2 else: - ww2 = energy_fourier(solver.comm, c.U_hat) / prod(params.N)**2 / 2 + ww2 = energy_fourier(solver.comm, c.U_hat) / prod(params.N) ** 2 / 2 kold[0] = kk if solver.rank == 0: @@ -127,11 +140,13 @@ def regression_test(context): U = solver.get_velocity(**context) curl = solver.get_curl(**context) w = solver.comm.reduce( - sum(curl.astype(float64) * curl.astype(float64)) / prod(params.N) / 2) + sum(curl.astype(float64) * curl.astype(float64)) / prod(params.N) / 2 + ) # Compute energy with double precision k = solver.comm.reduce( - sum(U.astype(float64) * U.astype(float64)) / prod(params.N) / 2) - config.solver.MemoryUsage('End', solver.comm) + sum(U.astype(float64) * U.astype(float64)) / prod(params.N) / 2 + ) + config.solver.MemoryUsage("End", solver.comm) if solver.rank == 0: assert round(asscalar(w) - 0.375249930801, params.ntol) == 0 assert round(asscalar(k) - 0.124953117517, params.ntol) == 0 @@ -139,35 +154,37 @@ def regression_test(context): if __name__ == "__main__": Re = 1e4 - U = 2 ** (1. / 3) - L = 1. + U = 2 ** (1.0 / 3) + L = 1.0 dt = 1e-12 config.update( { - 'nu': U * L / Re, # Viscosity - 'dt': dt, # Time step - 'T': 11*dt, # End time - 'L': [L, L, L], - 'M': [7, 7, 7], # Mesh size is pow(2, M[i]) in direction i + "nu": U * L / Re, # Viscosity + "dt": dt, # Time step + "T": 11 * dt, # End time + "L": [L, L, L], + "M": [7, 7, 7], # Mesh size is pow(2, M[i]) in direction i #'planner_effort': {'fft': 'FFTW_EXHAUSTIVE'}, #'decomposition': 'pencil', #'P1': 2 - }, "triplyperiodic" + }, + "triplyperiodic", ) config.triplyperiodic.add_argument("--compute_energy", type=int, default=2) config.triplyperiodic.add_argument("--plot_step", type=int, default=2) if plt is None: sol = get_solver(mesh="triplyperiodic") else: - sol = get_solver(update=update, regression_test=regression_test, - mesh="triplyperiodic") + sol = get_solver( + update=update, regression_test=regression_test, mesh="triplyperiodic" + ) context = sol.get_context() # Add curl to the stored results. For this we need to update the update_components # method used by the HDF5Writer class to compute the real fields that are stored WRITE_HDF5 = False - if config.params.solver == 'NS' and WRITE_HDF5: + if config.params.solver == "NS" and WRITE_HDF5: context.hdf5file.fname = "NS8.h5" context.hdf5file.components["curlx"] = context.curl[0] context.hdf5file.components["curly"] = context.curl[1] @@ -185,13 +202,14 @@ def update_components(**context): # Double check benchmark walltime start_time = time.time() - cProfile.runctx('solve(sol, context)', - globals(), locals(), 'profile.pstats') + cProfile.runctx("solve(sol, context)", globals(), locals(), "profile.pstats") end_time = time.time() - print('Run time: %f' % (end_time - start_time)) + print("Run time: %f" % (end_time - start_time)) - s = pstats.Stats('profile.pstats') - s.sort_stats('time').print_stats(12) + s = pstats.Stats("profile.pstats") + s.sort_stats("time").print_stats(12) - print('you can run:\n' - 'gprof2dot -f pstats profile.pstats | dot -Tpng -o profile.png') + print( + "you can run:\n" + "gprof2dot -f pstats profile.pstats | dot -Tpng -o profile.png" + ) diff --git a/bench/submit_benchmarks_legi.py b/bench/submit_benchmarks_legi.py index c10f177f..0db0a8a1 100644 --- a/bench/submit_benchmarks_legi.py +++ b/bench/submit_benchmarks_legi.py @@ -1,30 +1,33 @@ - - from fluiddyn.clusters.legi import Calcul7 as Cluster + # from fluiddyn.clusters.legi import Calcul8 as Cluster cluster = Cluster() -cluster.commands_setting_env.extend([ - 'module load p3dfft/2.7.4-mt', - 'module load pfft/1.0.6']) +cluster.commands_setting_env.extend( + ["module load p3dfft/2.7.4-mt", "module load pfft/1.0.6"] +) -def submit(nb_nodes, nb_cores_per_node=None, solver='ns2d'): +def submit(nb_nodes, nb_cores_per_node=None, solver="ns2d"): if nb_cores_per_node is None: nb_cores_per_node = cluster.nb_cores_per_node - nb_mpi = nb_cores_per_node*nb_nodes + nb_mpi = nb_cores_per_node * nb_nodes cluster.submit_command( - ('fluidsim bench 1024 -d 2 -s ' - '{} -o /.fsnet/data/legi/calcul9/home/augier3pi/fluidsim_bench' + ( + "fluidsim bench 1024 -d 2 -s " + "{} -o /.fsnet/data/legi/calcul9/home/augier3pi/fluidsim_bench" ).format(solver), - name_run=f'fluidsim-bench_{nb_mpi}', + name_run=f"fluidsim-bench_{nb_mpi}", nb_nodes=nb_nodes, # nb_cores_per_node=nb_cores_per_node, nb_cores_per_node=cluster.nb_cores_per_node, - walltime='00:04:00', - nb_mpi_processes=nb_mpi, omp_num_threads=1, + walltime="00:04:00", + nb_mpi_processes=nb_mpi, + omp_num_threads=1, ask=False, - delay_signal_walltime=None) + delay_signal_walltime=None, + ) + nb_nodes = 1 for nb_cores_per_node in [2, 4, 8, 12, 16, 20]: diff --git a/bench/util.py b/bench/util.py index 64f4bf45..e60ec421 100644 --- a/bench/util.py +++ b/bench/util.py @@ -6,85 +6,125 @@ import numpy as np -if os.getenv('SNIC_RESOURCE') is not None: +if os.getenv("SNIC_RESOURCE") is not None: # from fluiddyn.clusters.snic import ClusterSNIC as Cluster from fluiddyn.clusters.snic import Beskow36 as Cluster - cluster_type = 'snic' + + cluster_type = "snic" else: from fluiddyn.clusters.local import ClusterLocal as Cluster - cluster_type = 'local' + + cluster_type = "local" def create_common_params(n0, n1=None, n2=None): - params = Parameters('submit') - params._set_attrib('weak', False) - params._set_attrib('dry_run', False) - params._set_attrib('mode', '') - params._set_attrib('dim', 3) - params._set_attrib('shape', '') - params._set_attrib('output_dir', '') + params = Parameters("submit") + params._set_attrib("weak", False) + params._set_attrib("dry_run", False) + params._set_attrib("mode", "") + params._set_attrib("dim", 3) + params._set_attrib("shape", "") + params._set_attrib("output_dir", "") if n1 is None: n1 = n0 - params._set_child('two_d', dict( - shape=f'{n0} {n1}', time='00:20:00', - solver='ns2d', - fft_seq=['fft2d.with_fftw1d', - 'fft2d.with_fftw2d'], - fft=['fft2d.mpi_with_fftw1d', - 'fft2d.mpi_with_fftwmpi2d'], - nb_cores=np.array([2, 4, 8, 16, 32]), - nodes=[])) + params._set_child( + "two_d", + dict( + shape=f"{n0} {n1}", + time="00:20:00", + solver="ns2d", + fft_seq=["fft2d.with_fftw1d", "fft2d.with_fftw2d"], + fft=["fft2d.mpi_with_fftw1d", "fft2d.mpi_with_fftwmpi2d"], + nb_cores=np.array([2, 4, 8, 16, 32]), + nodes=[], + ), + ) if n2 is None: n2 = n0 - params._set_child('three_d', dict( - shape=f'{n0} {n1} {n2}', time='00:30:00', - solver='ns3d', - fft_seq=['fft3d.with_fftw3d'], - fft=['fft3d.mpi_with_fftw1d', - 'fft3d.mpi_with_fftwmpi3d', - 'fft3d.mpi_with_p3dfft', - 'fft3d.mpi_with_pfft'], - nb_cores=np.array([2, 4, 8, 16, 32]), - nodes=[])) + params._set_child( + "three_d", + dict( + shape=f"{n0} {n1} {n2}", + time="00:30:00", + solver="ns3d", + fft_seq=["fft3d.with_fftw3d"], + fft=[ + "fft3d.mpi_with_fftw1d", + "fft3d.mpi_with_fftwmpi3d", + "fft3d.mpi_with_p3dfft", + "fft3d.mpi_with_pfft", + ], + nb_cores=np.array([2, 4, 8, 16, 32]), + nodes=[], + ), + ) return params -def get_parser(prog='', description='', epilog=''): +def get_parser(prog="", description="", epilog=""): parser = argparse.ArgumentParser( prog=prog, description=description, formatter_class=argparse.RawDescriptionHelpFormatter, - epilog=epilog) - parser.add_argument('n0', nargs='?', type=int, default=None) - parser.add_argument('n1', nargs='?', type=int, default=None) - parser.add_argument('n2', nargs='?', type=int, default=None) + epilog=epilog, + ) + parser.add_argument("n0", nargs="?", type=int, default=None) + parser.add_argument("n1", nargs="?", type=int, default=None) + parser.add_argument("n2", nargs="?", type=int, default=None) parser.add_argument( - '-s', '--solver', type=str, default=None, - help='Any of the following solver keys: {}'.format( - available_solver_keys())) + "-s", + "--solver", + type=str, + default=None, + help="Any of the following solver keys: {}".format( + available_solver_keys() + ), + ) parser.add_argument( - '-d', '--dim', type=int, default=3, - help='dimension of the solver (default: 3)') + "-d", + "--dim", + type=int, + default=3, + help="dimension of the solver (default: 3)", + ) parser.add_argument( - '-n', '--dry-run', action='store_true', - help='simply print the commands which will be run') + "-n", + "--dry-run", + action="store_true", + help="simply print the commands which will be run", + ) parser.add_argument( - '-m', '--mode', default='seq-intra-inter', - help='could be "seq", "intra", "inter" or a combination of these') + "-m", + "--mode", + default="seq-intra-inter", + help='could be "seq", "intra", "inter" or a combination of these', + ) parser.add_argument( - '-nc', '--min-cores', type=int, default=1, - help='min. no. of processes to use (default: 1)') + "-nc", + "--min-cores", + type=int, + default=1, + help="min. no. of processes to use (default: 1)", + ) parser.add_argument( - '-nn', '--min-nodes', type=int, default=2, - help='max. no. of nodes to use for intra-node runs (default: 2)') + "-nn", + "--min-nodes", + type=int, + default=2, + help="max. no. of nodes to use for intra-node runs (default: 2)", + ) parser.add_argument( - '-xn', '--max-nodes', type=int, default=2, - help='max. no. of nodes to use for intra-node runs (default: 2)') + "-xn", + "--max-nodes", + type=int, + default=2, + help="max. no. of nodes to use for intra-node runs (default: 2)", + ) return parser @@ -103,61 +143,82 @@ def parser_to_params(parser): if args.min_cores > 1: log_min = np.log2(args.min_cores) params_dim.nb_cores = np.logspace( - log_min, 10, int(10 - log_min) + 1, base=2, dtype=int) + log_min, 10, int(10 - log_min) + 1, base=2, dtype=int + ) if args.max_nodes > 1: log_min = np.log2(args.min_nodes) log_max = np.log2(args.max_nodes) params_dim.nodes = np.logspace( - log_min, log_max, int(log_max - log_min) + 1, base=2, dtype=int) + log_min, log_max, int(log_max - log_min) + 1, base=2, dtype=int + ) params.dim = args.dim params.dry_run = args.dry_run params.mode = args.mode - params.shape = params_dim.shape.replace(' ', 'x') + params.shape = params_dim.shape.replace(" ", "x") return params, params_dim -def init_cluster(params, Cluster, prefix='snic', subdir='benchmarks'): +def init_cluster(params, Cluster, prefix="snic", subdir="benchmarks"): cluster = Cluster() - if cluster.name_cluster == 'beskow': - cluster.default_project = '2017-12-20' + if cluster.name_cluster == "beskow": + cluster.default_project = "2017-12-20" cluster.nb_cores_per_node = 32 else: - cluster.default_project = 'SNIC2017-12-20' + cluster.default_project = "SNIC2017-12-20" output_dir = params.output_dir = os.path.abspath( - './../../fluidsim-bench-results/{}/{}_{}'.format( - subdir, cluster.name_cluster, params.shape)) + "./../../fluidsim-bench-results/{}/{}_{}".format( + subdir, cluster.name_cluster, params.shape + ) + ) if not os.path.exists(output_dir): os.makedirs(output_dir) - print('Output directory: ', output_dir) - cluster.commands_unsetting_env.insert(0, 'fluidinfo -o ' + output_dir) + print("Output directory: ", output_dir) + cluster.commands_unsetting_env.insert(0, "fluidinfo -o " + output_dir) return cluster def submit( - params, params_dim, cluster, nb_nodes, nb_cores_per_node=None, - fft='all', cmd='bench'): + params, + params_dim, + cluster, + nb_nodes, + nb_cores_per_node=None, + fft="all", + cmd="bench", +): if nb_cores_per_node is None: nb_cores_per_node = cluster.nb_cores_per_node nb_mpi = nb_cores_per_node * nb_nodes nb_iterations = 20 - cmd = 'fluidsim {} -s {} {} -t {} -it {} -o {}'.format( - cmd, params_dim.solver, params_dim.shape, fft, nb_iterations, params.output_dir) + cmd = "fluidsim {} -s {} {} -t {} -it {} -o {}".format( + cmd, + params_dim.solver, + params_dim.shape, + fft, + nb_iterations, + params.output_dir, + ) if params.dry_run: - print('np =', nb_mpi, end=' ') + print("np =", nb_mpi, end=" ") print(cmd) else: cluster.submit_command( cmd, - name_run=f'{params_dim.solver}_{nb_mpi}', + name_run=f"{params_dim.solver}_{nb_mpi}", nb_nodes=nb_nodes, nb_cores_per_node=nb_cores_per_node, walltime=params_dim.time, - nb_mpi_processes=nb_mpi, omp_num_threads=1, - ask=False, bash=False, interactive=True, retain_script=True) + nb_mpi_processes=nb_mpi, + omp_num_threads=1, + ask=False, + bash=False, + interactive=True, + retain_script=True, + ) diff --git a/doc/example_from_numpy_doc.py b/doc/example_from_numpy_doc.py index 7009511b..731ea9b2 100644 --- a/doc/example_from_numpy_doc.py +++ b/doc/example_from_numpy_doc.py @@ -9,7 +9,7 @@ """ -import os # standard library imports first +import os # standard library imports first # Do NOT import using *, e.g. from numpy import * # @@ -34,7 +34,8 @@ from my_module import my_func, other_func -def foo(var1, var2, long_var_name='hi') : + +def foo(var1, var2, long_var_name="hi"): r"""A one-line summary that does not use variable names or the function name. diff --git a/doc/examples/internal_waves_coriolis/simul_idempotent.py b/doc/examples/internal_waves_coriolis/simul_idempotent.py index a6fd1e7e..93eb8a38 100644 --- a/doc/examples/internal_waves_coriolis/simul_idempotent.py +++ b/doc/examples/internal_waves_coriolis/simul_idempotent.py @@ -23,7 +23,13 @@ parser = argparse.ArgumentParser() parser.add_argument("amplitude", help="Amplitude of the movement", type=float) parser.add_argument("nz", help="Number of point over the z direction", type=int) -parser.add_argument("-me", "--max-elapsed", help="Maximum elapsed time", type=str, default="00:02:00") +parser.add_argument( + "-me", + "--max-elapsed", + help="Maximum elapsed time", + type=str, + default="00:02:00", +) args = parser.parse_args() mpi.printby0(args) @@ -53,7 +59,9 @@ if mpi.rank == 0: path_idempotent_file_exists = path_idempotent_file.exists() if mpi.nb_proc > 1: - path_idempotent_file_exists = mpi.comm.bcast(path_idempotent_file_exists, root=0) + path_idempotent_file_exists = mpi.comm.bcast( + path_idempotent_file_exists, root=0 + ) if not path_idempotent_file_exists: mpi.printby0("New simulation") diff --git a/doc/examples/internal_waves_coriolis/simul_ns3dstrat_waves.py b/doc/examples/internal_waves_coriolis/simul_ns3dstrat_waves.py index 8667c003..d1521cb8 100644 --- a/doc/examples/internal_waves_coriolis/simul_ns3dstrat_waves.py +++ b/doc/examples/internal_waves_coriolis/simul_ns3dstrat_waves.py @@ -15,9 +15,9 @@ # main input parameters N = 0.6 # rad/s -#if aspect_ratio == 4: +# if aspect_ratio == 4: omega_f = 0.807 * N # rad/s -#elif aspect_ratio == 6: +# elif aspect_ratio == 6: # omega_f = 0.674 * N # rad/s delta_omega_f = 0.1 * omega_f # rad/s @@ -133,4 +133,4 @@ ) -sim.output.spect_energy_budg.plot_fluxes(key_k='kh') +sim.output.spect_energy_budg.plot_fluxes(key_k="kh") diff --git a/doc/examples/internal_waves_coriolis/simul_ns3dstrat_waves_forcinginscript.py b/doc/examples/internal_waves_coriolis/simul_ns3dstrat_waves_forcinginscript.py index a7c76d6c..517ae311 100644 --- a/doc/examples/internal_waves_coriolis/simul_ns3dstrat_waves_forcinginscript.py +++ b/doc/examples/internal_waves_coriolis/simul_ns3dstrat_waves_forcinginscript.py @@ -174,9 +174,7 @@ def step_func(x): def compute_forcing_each_time(self): - """This function is called by the forcing_maker to compute the forcing - - """ + """This function is called by the forcing_maker to compute the forcing""" sim = self.sim time = sim.time_stepping.t % period_forcing coef_forcing_time_x = calcul_forcing_time_x(time) diff --git a/doc/examples/internal_waves_coriolis/submit_simul.py b/doc/examples/internal_waves_coriolis/submit_simul.py index 4df5edb8..64bf4ca4 100644 --- a/doc/examples/internal_waves_coriolis/submit_simul.py +++ b/doc/examples/internal_waves_coriolis/submit_simul.py @@ -3,14 +3,14 @@ cluster = Cluster() cluster.commands_setting_env = [ - 'source /etc/profile', + "source /etc/profile", 'export PATH="$HOME/miniconda3/bin:$PATH"' - 'export FLUIDSIM_PATH=/fsnet/project/watu/2019/19INTSIM/sim_data' + "export FLUIDSIM_PATH=/fsnet/project/watu/2019/19INTSIM/sim_data", ] cluster.submit_script( - 'simul_ns3dstrat_waves.py', - name_run='fld_igw3d', + "simul_ns3dstrat_waves.py", + name_run="fld_igw3d", nb_cores_per_node=10, nb_mpi_processes=10, omp_num_threads=1, diff --git a/doc/examples/launch_simul.py b/doc/examples/launch_simul.py index b052e9d7..a270cc8f 100644 --- a/doc/examples/launch_simul.py +++ b/doc/examples/launch_simul.py @@ -1,5 +1,5 @@ - from fluiddyn.clusters.legi import Calcul as Cluster + # or # from fluiddyn.clusters.legi import Calcul7 as Cluster # or @@ -8,13 +8,13 @@ cluster = Cluster() cluster.commands_setting_env = [ - 'source /etc/profile', - 'export PATH="$HOME/miniconda3/bin:$PATH"' + "source /etc/profile", + 'export PATH="$HOME/miniconda3/bin:$PATH"', ] cluster.submit_script( - 'simul_ns2d.py', - name_run='fld_example', + "simul_ns2d.py", + name_run="fld_example", nb_cores_per_node=4, nb_mpi_processes=4, omp_num_threads=1, diff --git a/doc/examples/milestone/run_simul.py b/doc/examples/milestone/run_simul.py index c3d48dc0..d4da9e71 100644 --- a/doc/examples/milestone/run_simul.py +++ b/doc/examples/milestone/run_simul.py @@ -50,7 +50,11 @@ ) parser.add_argument( - "-np", "--n_periods", type=int, default=5, help="Number of periods", + "-np", + "--n_periods", + type=int, + default=5, + help="Number of periods", ) parser.add_argument( diff --git a/doc/examples/simul_ad1d.py b/doc/examples/simul_ad1d.py index 4ca93e61..5f9a66ba 100644 --- a/doc/examples/simul_ad1d.py +++ b/doc/examples/simul_ad1d.py @@ -1,31 +1,30 @@ - from fluidsim.solvers.ad1d.solver import Simul params = Simul.create_default_params() -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" -params.U = 1. +params.U = 1.0 params.oper.nx = 200 -params.oper.Lx = 1. +params.oper.Lx = 1.0 -params.time_stepping.type_time_scheme = 'RK2' +params.time_stepping.type_time_scheme = "RK2" params.nu_2 = 0.01 params.time_stepping.t_end = 0.4 params.time_stepping.USE_CFL = True -params.init_fields.type = 'gaussian' +params.init_fields.type = "gaussian" params.output.periods_print.print_stdout = 0.25 params.output.periods_save.phys_fields = 0.1 -params.output.periods_plot.phys_fields = 0. +params.output.periods_plot.phys_fields = 0.0 -params.output.phys_fields.field_to_plot = 's' +params.output.phys_fields.field_to_plot = "s" sim = Simul(params) diff --git a/doc/examples/simul_ns2d.py b/doc/examples/simul_ns2d.py index 900f4312..9bdb4661 100644 --- a/doc/examples/simul_ns2d.py +++ b/doc/examples/simul_ns2d.py @@ -8,20 +8,20 @@ params.oper.Lx = params.oper.Ly = Lh = 2 * pi delta_x = Lh / nh -params.nu_8 = 2.*params.forcing.forcing_rate**(1./3)*delta_x**8 +params.nu_8 = 2.0 * params.forcing.forcing_rate ** (1.0 / 3) * delta_x ** 8 -params.time_stepping.t_end = 2. +params.time_stepping.t_end = 2.0 -params.init_fields.type = 'dipole' +params.init_fields.type = "dipole" params.forcing.enable = True -params.forcing.type = 'proportional' +params.forcing.type = "proportional" -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" params.output.periods_print.print_stdout = 0.25 -params.output.periods_save.phys_fields = 1. +params.output.periods_save.phys_fields = 1.0 params.output.periods_save.spectra = 0.5 params.output.periods_save.spatial_means = 0.05 params.output.periods_save.spect_energy_budg = 0.5 diff --git a/doc/examples/simul_ns2d_plot.py b/doc/examples/simul_ns2d_plot.py index f46577aa..eba020ad 100644 --- a/doc/examples/simul_ns2d_plot.py +++ b/doc/examples/simul_ns2d_plot.py @@ -1,29 +1,28 @@ - import fluiddyn as fld from fluidsim.solvers.ns2d.solver import Simul params = Simul.create_default_params() -params.short_name_type_run = 'test' +params.short_name_type_run = "test" params.oper.nx = params.oper.ny = nh = 32 params.oper.Lx = params.oper.Ly = Lh = 10 delta_x = Lh / nh -params.nu_8 = 2.*params.forcing.forcing_rate**(1./3)*delta_x**8 +params.nu_8 = 2.0 * params.forcing.forcing_rate ** (1.0 / 3) * delta_x ** 8 -params.time_stepping.t_end = 10. +params.time_stepping.t_end = 10.0 -params.init_fields.type = 'dipole' +params.init_fields.type = "dipole" params.forcing.enable = True -params.forcing.type = 'tcrandom' +params.forcing.type = "tcrandom" -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" params.output.periods_print.print_stdout = 0.5 -params.output.periods_save.phys_fields = 1. +params.output.periods_save.phys_fields = 1.0 params.output.periods_save.spectra = 0.5 params.output.periods_save.spatial_means = 0.05 params.output.periods_save.spect_energy_budg = 0.5 @@ -36,10 +35,10 @@ params.output.spect_energy_budg.HAS_TO_PLOT_SAVED = True params.output.increments.HAS_TO_PLOT_SAVED = True -params.output.phys_fields.field_to_plot = 'rot' +params.output.phys_fields.field_to_plot = "rot" sim = Simul(params) - + sim.time_stepping.start() sim.output.phys_fields.plot() diff --git a/doc/examples/simul_ns2d_plotphys.py b/doc/examples/simul_ns2d_plotphys.py index f368a7e0..b7e3a448 100644 --- a/doc/examples/simul_ns2d_plotphys.py +++ b/doc/examples/simul_ns2d_plotphys.py @@ -1,25 +1,24 @@ - import fluiddyn as fld from fluidsim.solvers.ns2d.solver import Simul params = Simul.create_default_params() -params.short_name_type_run = 'test' +params.short_name_type_run = "test" -params.oper.nx = params.oper.ny = nh = 32*2 -params.oper.Lx = params.oper.Ly = Lh = 10. +params.oper.nx = params.oper.ny = nh = 32 * 2 +params.oper.Lx = params.oper.Ly = Lh = 10.0 delta_x = Lh / nh -params.nu_8 = 2e-3*params.forcing.forcing_rate**(1./3)*delta_x**8 +params.nu_8 = 2e-3 * params.forcing.forcing_rate ** (1.0 / 3) * delta_x ** 8 -params.time_stepping.t_end = 10. +params.time_stepping.t_end = 10.0 -params.init_fields.type = 'dipole' +params.init_fields.type = "dipole" params.forcing.enable = True -params.forcing.type = 'tcrandom' +params.forcing.type = "tcrandom" -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" params.output.periods_plot.phys_fields = 0.1 params.output.periods_save.phys_fields = 0.2 @@ -31,10 +30,12 @@ sim = Simul(params) sim.time_stepping.start() -print(""" +print( + """ A movie can be produced with the command (using ffmpeg): sim.output.phys_fields.animate(dt_frame_in_sec=0.1, dt_equations=0.08, repeat=False, save_file=1, tmax=10) -""") +""" +) fld.show() diff --git a/doc/examples/simul_ns2d_restart.py b/doc/examples/simul_ns2d_restart.py index a642b129..bb0e9937 100644 --- a/doc/examples/simul_ns2d_restart.py +++ b/doc/examples/simul_ns2d_restart.py @@ -1,4 +1,3 @@ - from fluidsim import load_for_restart, path_dir_results path_dir_root = path_dir_results / "examples" diff --git a/doc/examples/simul_ns2dbouss_initfields_in_script.py b/doc/examples/simul_ns2dbouss_initfields_in_script.py index 4439097c..a4239fdf 100644 --- a/doc/examples/simul_ns2dbouss_initfields_in_script.py +++ b/doc/examples/simul_ns2dbouss_initfields_in_script.py @@ -13,21 +13,21 @@ params = Simul.create_default_params() -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" params.oper.nx = nx = 64 -params.oper.ny = nx//2 +params.oper.ny = nx // 2 params.oper.Lx = lx = 2 params.oper.Ly = lx params.oper.coef_dealiasing = 0.7 params.nu_8 = 1e-10 -params.time_stepping.t_end = 5. +params.time_stepping.t_end = 5.0 -params.init_fields.type = 'in_script' +params.init_fields.type = "in_script" -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" params.output.periods_print.print_stdout = 0.5 params.output.periods_save.phys_fields = 0.1 params.output.periods_save.spatial_means = 0.1 @@ -38,10 +38,10 @@ rot = 1e-6 * sim.oper.create_arrayX_random() X = sim.oper.X Y = sim.oper.Y -x0 = y0 = 1. -R2 = (X-x0)**2 + (Y-y0)**2 +x0 = y0 = 1.0 +R2 = (X - x0) ** 2 + (Y - y0) ** 2 r0 = 0.2 -b = -np.exp(-R2/r0**2) +b = -np.exp(-R2 / r0 ** 2) sim.state.init_from_rotb(rot, b) # In this case (params.init_fields.type = 'in_script') if we want to plot the @@ -55,9 +55,9 @@ if rank == 0: print( - '\nTo display a video of this simulation, you can do:\n' - f'cd {sim.output.path_run}' + - """ + "\nTo display a video of this simulation, you can do:\n" + f"cd {sim.output.path_run}" + + """ ipython # then in ipython (copy the 3 lines in the terminal): @@ -66,4 +66,5 @@ sim = load_sim_for_plot() sim.output.phys_fields.animate('b', dt_frame_in_sec=0.3, dt_equations=0.1) -""") +""" + ) diff --git a/doc/examples/simul_ns2dbouss_staticlayer.py b/doc/examples/simul_ns2dbouss_staticlayer.py index 1be85139..e25f01ed 100644 --- a/doc/examples/simul_ns2dbouss_staticlayer.py +++ b/doc/examples/simul_ns2dbouss_staticlayer.py @@ -18,11 +18,11 @@ params = Simul.create_default_params() -params.output.sub_directory = 'examples' -params.short_name_type_run = 'staticlayer' +params.output.sub_directory = "examples" +params.short_name_type_run = "staticlayer" params.oper.nx = nx = 128 -params.oper.ny = ny = nx//4 +params.oper.ny = ny = nx // 4 params.oper.Lx = lx = 4 params.oper.Ly = ly = lx * ny / nx @@ -31,18 +31,18 @@ params.nu_8 = 1e-10 -params.time_stepping.t_end = 8. +params.time_stepping.t_end = 8.0 # we need small time step for a strong forcing params.time_stepping.USE_CFL = False params.time_stepping.deltat0 = 0.02 -params.init_fields.type = 'in_script' +params.init_fields.type = "in_script" params.forcing.enable = True -params.forcing.type = 'in_script_coarse' +params.forcing.type = "in_script_coarse" params.forcing.nkmax_forcing = 10 -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" params.output.periods_print.print_stdout = 0.5 params.output.periods_save.phys_fields = 0.1 params.output.periods_save.spatial_means = 0.1 @@ -59,7 +59,7 @@ y0 = ly / 2 blob_height = 0.4 blob_width = 0.5 -b = -np.exp(-(X - x0)**2/blob_width**2 - (Y - y0)**2/blob_height**2) +b = -np.exp(-((X - x0) ** 2) / blob_width ** 2 - (Y - y0) ** 2 / blob_height ** 2) sim.state.init_from_rotb(rot, b) # in this case (params.init_fields.type = 'manual') if we want to plot the @@ -75,8 +75,8 @@ forcing_maker = sim.forcing.forcing_maker oper = forcing_maker.oper_coarse Y = oper.Y - d = ly/6 - alpha = - 80 * (np.exp(-Y**2/d**2) + np.exp(-(Y - ly)**2/d**2)) + d = ly / 6 + alpha = -80 * (np.exp(-(Y ** 2) / d ** 2) + np.exp(-((Y - ly) ** 2) / d ** 2)) # on-the-fly plot has_to_animate = 1 @@ -84,12 +84,12 @@ # initialization of the on-the-fly plot fig = plt.figure(figsize=(10, 10)) subplots = fig.subplots(2, 3) - title = fig.suptitle('time = 0') + title = fig.suptitle("time = 0") xs = oper.x ys = oper.y pmeshs = np.empty_like(subplots) - subtitles = (('ux', 'uy', 'rot'), ('fx', 'fy', 'frot')) + subtitles = (("ux", "uy", "rot"), ("fx", "fy", "frot")) for (i0, i1), ax in np.ndenumerate(subplots): pmeshs[i0, i1] = ax.pcolormesh(xs, ys, alpha) @@ -102,12 +102,11 @@ # sys.exit() def compute_forcingc_fft_each_time(self): - """This function is called by the forcing_maker to compute the forcing - - """ - rot_fft = self.sim.state.state_spect.get_var('rot_fft') + """This function is called by the forcing_maker to compute the forcing""" + rot_fft = self.sim.state.state_spect.get_var("rot_fft") rot_fft = self.oper.coarse_seq_from_fft_loc( - rot_fft, self.shapeK_loc_coarse) + rot_fft, self.shapeK_loc_coarse + ) ux_fft, uy_fft = oper.vecfft_from_rotfft(rot_fft) ux = oper.ifft(ux_fft) uy = oper.ifft(uy_fft) @@ -131,7 +130,7 @@ def compute_forcingc_fft_each_time(self): pmesh.set_array(arr[:-1, :-1].ravel()) pmesh.set_clim(arr.min(), arr.max()) - title.set_text(f'time {self.sim.time_stepping.t:.2f}') + title.set_text(f"time {self.sim.time_stepping.t:.2f}") fig.canvas.draw() plt.pause(1e-4) @@ -139,16 +138,17 @@ def compute_forcingc_fft_each_time(self): return frot_fft forcing_maker.monkeypatch_compute_forcingc_fft_each_time( - compute_forcingc_fft_each_time) + compute_forcingc_fft_each_time + ) # and finally time stepping sim.time_stepping.start() if rank == 0: print( - '\nTo display a video of this simulation, you can do:\n' - f'cd {sim.output.path_run}' + - """ + "\nTo display a video of this simulation, you can do:\n" + f"cd {sim.output.path_run}" + + """ ipython # then in ipython (copy the 3 lines in the terminal): @@ -157,6 +157,7 @@ def compute_forcingc_fft_each_time(self): sim = load_sim_for_plot() sim.output.phys_fields.animate('uy', dt_frame_in_sec=0.3, dt_equations=0.1) -""") +""" + ) plt.show() diff --git a/doc/examples/simul_ns2dstrat_forcing_in_script.py b/doc/examples/simul_ns2dstrat_forcing_in_script.py index f5637d60..a3d7f252 100644 --- a/doc/examples/simul_ns2dstrat_forcing_in_script.py +++ b/doc/examples/simul_ns2dstrat_forcing_in_script.py @@ -13,36 +13,36 @@ params = Simul.create_default_params() -params.output.sub_directory = 'examples' -params.short_name_type_run = 'forcinginscript' +params.output.sub_directory = "examples" +params.short_name_type_run = "forcinginscript" params.oper.nx = nx = 64 -params.oper.ny = nx//2 +params.oper.ny = nx // 2 params.oper.Lx = lx = 2 * pi -params.oper.Ly = lx/2 +params.oper.Ly = lx / 2 params.oper.coef_dealiasing = 0.7 params.nu_8 = 1e-9 -params.N = 1. # Brunt Vaisala frequency +params.N = 1.0 # Brunt Vaisala frequency -params.time_stepping.t_end = 10. +params.time_stepping.t_end = 10.0 -params.init_fields.type = 'noise' +params.init_fields.type = "noise" params.forcing.enable = True -params.forcing.type = 'in_script_coarse' +params.forcing.type = "in_script_coarse" params.forcing.nkmax_forcing = 12 -params.forcing.key_forced = 'rot_fft' +params.forcing.key_forced = "rot_fft" -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" params.output.periods_print.print_stdout = 0.5 params.output.periods_save.phys_fields = 0.2 params.output.periods_save.spectra = 0.5 params.output.periods_save.spatial_means = 0.05 -params.output.periods_save.spect_energy_budg = 1. -params.output.periods_save.increments = 1. +params.output.periods_save.spect_energy_budg = 1.0 +params.output.periods_save.increments = 1.0 sim = Simul(params) @@ -50,23 +50,24 @@ if rank == 0: forcing_maker = sim.forcing.forcing_maker oper = forcing_maker.oper_coarse - forcing0 = 2*np.cos(2*pi*oper.Y / oper.ly) - omega = 2*pi + forcing0 = 2 * np.cos(2 * pi * oper.Y / oper.ly) + omega = 2 * pi def compute_forcingc_each_time(self): - return forcing0 * np.sin(omega*sim.time_stepping.t) + return forcing0 * np.sin(omega * sim.time_stepping.t) forcing_maker.monkeypatch_compute_forcingc_each_time( - compute_forcingc_each_time) + compute_forcingc_each_time + ) sim.time_stepping.start() if rank == 0: print( - '\nTo display a video of this simulation, you can do:\n' - f'cd {sim.output.path_run}' + - """ + "\nTo display a video of this simulation, you can do:\n" + f"cd {sim.output.path_run}" + + """ ipython # then in ipython (copy the 3 lines in the terminal): @@ -75,4 +76,5 @@ def compute_forcingc_each_time(self): sim = load_sim_for_plot() sim.output.phys_fields.animate('b', dt_frame_in_sec=0.1, dt_equations=0.1) -""") +""" + ) diff --git a/doc/examples/simul_ns3dbouss_staticlayer.py b/doc/examples/simul_ns3dbouss_staticlayer.py index b68b4aaa..9169d121 100644 --- a/doc/examples/simul_ns3dbouss_staticlayer.py +++ b/doc/examples/simul_ns3dbouss_staticlayer.py @@ -16,21 +16,21 @@ params = Simul.create_default_params() -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" -nx = 144//2 +nx = 144 // 2 ny = nx -nz = nx//2 +nz = nx // 2 lz = 2 params.oper.nx = nx params.oper.ny = ny params.oper.nz = nz -params.oper.Lx = lx = lz/nz*nx -params.oper.Ly = ly = lz/nz*ny +params.oper.Lx = lx = lz / nz * nx +params.oper.Ly = ly = lz / nz * ny params.oper.Lz = lz # rotation -params.f = 1. +params.f = 1.0 r""" @@ -56,23 +56,23 @@ """ n = 8 -C = 1. -dx = lx/nx +C = 1.0 +dx = lx / nx B = 1 D = 1 -eps = 1e-2*B**(3/2)*D**(1/2) -params.nu_8 = (dx/C)**((3*n-2)/3) * eps**(1/3) +eps = 1e-2 * B ** (3 / 2) * D ** (1 / 2) +params.nu_8 = (dx / C) ** ((3 * n - 2) / 3) * eps ** (1 / 3) # printby0(f'nu_8 = {params.nu_8:.3e}') params.time_stepping.USE_T_END = True -params.time_stepping.t_end = 4. +params.time_stepping.t_end = 4.0 params.time_stepping.deltat_max = 0.01 -params.init_fields.type = 'in_script' +params.init_fields.type = "in_script" params.forcing.enable = True -params.forcing.type = 'in_script' +params.forcing.type = "in_script" params.output.periods_print.print_stdout = 1e-1 @@ -82,28 +82,34 @@ sim = Simul(params) # here we have to initialize the flow fields -variables = {k: 1e-2 * sim.oper.create_arrayX_random() - for k in ('vx', 'vy', 'vz')} +variables = { + k: 1e-2 * sim.oper.create_arrayX_random() for k in ("vx", "vy", "vz") +} X, Y, Z = sim.oper.get_XYZ_loc() -width_step = max(4*dx, 0.1) +width_step = max(4 * dx, 0.1) + def step_func(x): """Activation function""" - return 0.5*(np.tanh(x/width_step) + 1) + return 0.5 * (np.tanh(x / width_step) + 1) -x0 = lx/2. -y0 = ly/2. -Rh2 = (X-x0)**2 + (Y-y0)**2 +x0 = lx / 2.0 +y0 = ly / 2.0 +Rh2 = (X - x0) ** 2 + (Y - y0) ** 2 r0 = 0.5 -d_forcing = lz/10 +d_forcing = lz / 10 d_forcing_b = 1.2 * d_forcing -b = -0.5*(1-np.tanh((Rh2 - r0**2)/0.2**2)) * \ - step_func(Z - d_forcing_b) * step_func(-(Z - (lz - d_forcing_b))) -variables['b'] = b +b = ( + -0.5 + * (1 - np.tanh((Rh2 - r0 ** 2) / 0.2 ** 2)) + * step_func(Z - d_forcing_b) + * step_func(-(Z - (lz - d_forcing_b))) +) +variables["b"] = b sim.state.init_statephys_from(**variables) @@ -127,31 +133,32 @@ def step_func(x): """ gamma = 2 n_dt = 4 -coef_sigma = gamma/n_dt +coef_sigma = gamma / n_dt def compute_forcing_fft_each_time(self): - """This function is called by the forcing_maker to compute the forcing + """This function is called by the forcing_maker to compute the forcing""" - """ def compute_forcing_1var(key): vi = self.sim.state.state_phys.get_var(key) - sigma = coef_sigma/self.sim.time_stepping.deltat - fi = - sigma * alpha * vi + sigma = coef_sigma / self.sim.time_stepping.deltat + fi = -sigma * alpha * vi return oper.fft(fi) - keys = ('vz',) - result = {key + '_fft': compute_forcing_1var(key) - for key in keys} + keys = ("vz",) + result = {key + "_fft": compute_forcing_1var(key) for key in keys} return result + sim.forcing.forcing_maker.monkeypatch_compute_forcing_fft_each_time( - compute_forcing_fft_each_time) + compute_forcing_fft_each_time +) # finally we start the simulation sim.time_stepping.start() -printby0(f""" +printby0( + f""" # To visualize the output with Paraview, create a file states_phys.xmf with: fluidsim-create-xml-description {sim.output.path_run} @@ -168,4 +175,5 @@ def compute_forcing_1var(key): sim.output.phys_fields.set_equation_crosssection('x={lx/2}') sim.output.phys_fields.animate('b') -""") +""" +) diff --git a/doc/examples/simul_ns3dstrat_check_energy_conserv.py b/doc/examples/simul_ns3dstrat_check_energy_conserv.py index a9c2611a..063d54d7 100644 --- a/doc/examples/simul_ns3dstrat_check_energy_conserv.py +++ b/doc/examples/simul_ns3dstrat_check_energy_conserv.py @@ -12,9 +12,9 @@ params = Simul.create_default_params() -params.f = 1. +params.f = 1.0 -params.output.sub_directory = 'examples' +params.output.sub_directory = "examples" nx = 48 ny = 64 @@ -24,10 +24,10 @@ params.oper.ny = ny params.oper.nz = nz params.oper.Lx = Lx -params.oper.Ly = Ly = Lx/nx*ny -params.oper.Lz = Lz = Lx/nx*nz +params.oper.Ly = Ly = Lx / nx * ny +params.oper.Lz = Lz = Lx / nx * nz -params.short_name_type_run = 'checkenergy' +params.short_name_type_run = "checkenergy" params.oper.coef_dealiasing = 0.5 @@ -36,10 +36,10 @@ params.time_stepping.it_end = 10 params.time_stepping.deltat0 = 1e-2 -params.init_fields.type = 'noise' +params.init_fields.type = "noise" params.output.periods_print.print_stdout = 1e-10 -params.output.periods_save.phys_fields = 0. +params.output.periods_save.phys_fields = 0.0 sim = Simul(params) diff --git a/doc/examples/simul_ns3dstrat_forcing_random.py b/doc/examples/simul_ns3dstrat_forcing_random.py index 111fba12..a50d85d5 100644 --- a/doc/examples/simul_ns3dstrat_forcing_random.py +++ b/doc/examples/simul_ns3dstrat_forcing_random.py @@ -1,4 +1,3 @@ - from fluiddyn.util.mpi import printby0 from fluidsim.solvers.ns3d.strat.solver import Simul @@ -30,7 +29,7 @@ printby0(f"nu_8 = {params.nu_8:.3e}") params.init_fields.type = "noise" -params.init_fields.noise.length = 1. +params.init_fields.noise.length = 1.0 params.init_fields.noise.velo_max = 0.1 params.forcing.enable = True diff --git a/doc/examples/time_stepping/3d/run_simul.py b/doc/examples/time_stepping/3d/run_simul.py index f2f72943..7ce13093 100644 --- a/doc/examples/time_stepping/3d/run_simul.py +++ b/doc/examples/time_stepping/3d/run_simul.py @@ -66,7 +66,10 @@ ) parser.add_argument( - "--t_end", help="params.time_stepping.t_end", type=float, default=20.0, + "--t_end", + help="params.time_stepping.t_end", + type=float, + default=20.0, ) parser.add_argument( @@ -84,7 +87,10 @@ ) parser.add_argument( - "--Re", help="Reynolds number", type=float, default=1600.0, + "--Re", + help="Reynolds number", + type=float, + default=1600.0, ) @@ -92,11 +98,17 @@ L = 1 parser.add_argument( - "--velo_max_noise", help="noise level", type=float, default=0.0, + "--velo_max_noise", + help="noise level", + type=float, + default=0.0, ) parser.add_argument( - "--length_noise", help="noise length", type=float, default=L / 5, + "--length_noise", + help="noise length", + type=float, + default=L / 5, ) parser.add_argument( diff --git a/scripts/bugs/no_energy_conservation_ns2dstrat.py b/scripts/bugs/no_energy_conservation_ns2dstrat.py index c5989d03..8b229dcc 100644 --- a/scripts/bugs/no_energy_conservation_ns2dstrat.py +++ b/scripts/bugs/no_energy_conservation_ns2dstrat.py @@ -36,7 +36,7 @@ params.oper.nx = nx = 64 params.oper.ny = ny = nx -params.oper.Lx = 10. +params.oper.Lx = 10.0 params.oper.Ly = params.oper.Lx * (ny / nx) # params.oper.coef_dealiasing = 0.5 @@ -50,17 +50,17 @@ # params.nu_8 = 1e-14 # params.N = 3. -params.init_fields.type = 'noise' +params.init_fields.type = "noise" params.init_fields.noise.length = 1 params.init_fields.noise.velo_max = 1e-4 params.forcing.enable = True # params.forcing.type = 'tcrandom_anisotropic' -params.forcing.type = 'tcrandom' +params.forcing.type = "tcrandom" params.forcing.nkmax_forcing = 10 params.forcing.nkmin_forcing = 4 -params.forcing.tcrandom.time_correlation = 1. +params.forcing.tcrandom.time_correlation = 1.0 # params.forcing.key_forced = 'ap_fft' # params.forcing.normalized.which_root = 'minabs' diff --git a/scripts/bugs/wrong_output_it0_initfield.py b/scripts/bugs/wrong_output_it0_initfield.py index 70dc17d5..2041c717 100644 --- a/scripts/bugs/wrong_output_it0_initfield.py +++ b/scripts/bugs/wrong_output_it0_initfield.py @@ -29,18 +29,18 @@ params.oper.nx = nx = 128 params.oper.ny = ny = nx // 2 -params.oper.Lx = 10. +params.oper.Lx = 10.0 params.oper.Ly = params.oper.Lx * (ny / nx) params.nu_8 = 1e-10 -params.init_fields.type = 'noise' -params.init_fields.noise.velo_max = 1. -params.init_fields.noise.length = 1. +params.init_fields.type = "noise" +params.init_fields.noise.velo_max = 1.0 +params.init_fields.noise.length = 1.0 params.init_fields.modif_after_init = True params.forcing.enable = True -params.forcing.type = 'tcrandom' +params.forcing.type = "tcrandom" params.forcing.nkmax_forcing = 10 params.forcing.nkmin_forcing = 4 @@ -57,7 +57,7 @@ sim = Simul(params) # we modify the state after the init... -sim.state.state_phys *= 2. +sim.state.state_phys *= 2.0 sim.state.statespect_from_statephys() sim.time_stepping.start() diff --git a/scripts/launch/simul_test_conserveE.py b/scripts/launch/simul_test_conserveE.py index 825a270b..e0db187a 100644 --- a/scripts/launch/simul_test_conserveE.py +++ b/scripts/launch/simul_test_conserveE.py @@ -1,10 +1,10 @@ # !/usr/bin/env python -#coding=utf8 +# coding=utf8 # # run simul_test_conserveE.py # mpirun -np 8 python simul_test_conserveE.py -key_solver = 'NS2D' +key_solver = "NS2D" # key_solver = 'SW1l' # key_solver = 'SW1l.onlywaves' # key_solver = 'SW1l.exactlin' @@ -16,25 +16,25 @@ solver = fluidsim.import_module_solver_from_key(key_solver) params = solver.Simul.create_default_params() -params.short_name_type_run = 'conservE' +params.short_name_type_run = "conservE" -nh = 16*4*8 -Lh = 2*np.pi +nh = 16 * 4 * 8 +Lh = 2 * np.pi params.oper.nx = nh params.oper.ny = nh params.oper.Lx = Lh params.oper.Ly = Lh -params.oper.coef_dealiasing = 2./3 +params.oper.coef_dealiasing = 2.0 / 3 -params.nu_8 = 0. -params.nu_4 = 0. -params.nu_2 = 0. -params.nu_m4 = 0. +params.nu_8 = 0.0 +params.nu_4 = 0.0 +params.nu_2 = 0.0 +params.nu_m4 = 0.0 try: - params.f = 1. - params.c2 = 200. + params.f = 1.0 + params.c2 = 200.0 except KeyError: pass @@ -43,32 +43,32 @@ params.time_stepping.it_end = 5 params.time_stepping.USE_T_END = False -params.time_stepping.type_time_scheme = 'RK4' +params.time_stepping.type_time_scheme = "RK4" -params.init_fields.type_flow_init = 'DIPOLE' +params.init_fields.type_flow_init = "DIPOLE" -params.output.periods_print.print_stdout = 10.e-15 +params.output.periods_print.print_stdout = 10.0e-15 -params.output.periods_save.phys_fields = 0. -params.output.periods_save.spectra = 0. -params.output.periods_save.spect_energy_budg = 0. -params.output.periods_save.increments = 0. +params.output.periods_save.phys_fields = 0.0 +params.output.periods_save.spectra = 0.0 +params.output.periods_save.spect_energy_budg = 0.0 +params.output.periods_save.increments = 0.0 try: - params.output.periods_save.pdf = 0. + params.output.periods_save.pdf = 0.0 params.output.periods_save.time_signals_fft = False except KeyError: pass -params.output.periods_plot.phys_fields = 0. +params.output.periods_plot.phys_fields = 0.0 -params.time_stepping.deltat0 = 1.e-1 +params.time_stepping.deltat0 = 1.0e-1 sim = solver.Simul(params) sim.time_stepping.start() -params.time_stepping.deltat0 = 1.e-2 +params.time_stepping.deltat0 = 1.0e-2 sim = solver.Simul(params) sim.time_stepping.start() diff --git a/scripts/launch/sw1l_force_vortex_grid.py b/scripts/launch/sw1l_force_vortex_grid.py index 19a6a493..473da689 100755 --- a/scripts/launch/sw1l_force_vortex_grid.py +++ b/scripts/launch/sw1l_force_vortex_grid.py @@ -6,8 +6,8 @@ from fluidsim.util.util import mpi params = Simul.create_default_params() -params.short_name_type_run = 'vortex_grid' -params.output.sub_directory = 'beskow_tests' +params.short_name_type_run = "vortex_grid" +params.output.sub_directory = "beskow_tests" # --------Grid parameters--------- params.oper.nx = params.oper.ny = nh = 512 @@ -15,33 +15,33 @@ delta_x = Lh / nh # -----Numerical method parameters---- -params.f = 0. +params.f = 0.0 params.c2 = 100 -params.oper.coef_dealiasing = 2./3 -params.time_stepping.t_end = 100. -params.init_fields.type = 'vortex_grid' +params.oper.coef_dealiasing = 2.0 / 3 +params.time_stepping.t_end = 100.0 +params.init_fields.type = "vortex_grid" # ----- Forcing parameters ----------- params.FORCING = True -params.forcing.type = 'waves' +params.forcing.type = "waves" params.forcing.nkmax_forcing = 8 params.forcing.nkmin_forcing = 5 -params.forcing.forcing_rate = 1. +params.forcing.forcing_rate = 1.0 -#k_max = pi / delta_x * params.oper.coef_dealiasing # Smallest resolved scale -#k_d = k_max / 2.5 -#length_scale = pi / k_d -#params.nu_8 = params.forcing.forcing_rate**(1./3) * length_scale ** (22./3) +# k_max = pi / delta_x * params.oper.coef_dealiasing # Smallest resolved scale +# k_d = k_max / 2.5 +# length_scale = pi / k_d +# params.nu_8 = params.forcing.forcing_rate**(1./3) * length_scale ** (22./3) # ------Save file / StdOut parameters------------ params.output.periods_print.print_stdout = 0.1 params.output.period_refresh_plots = 0.5 params.output.periods_plot.phys_fields = 0.0 -#params.output.spatial_means.HAS_TO_PLOT_SAVED = True -#params.output.spectra.HAS_TO_PLOT_SAVED = True -#params.output.spect_energy_budg.HAS_TO_PLOT_SAVED = True +# params.output.spatial_means.HAS_TO_PLOT_SAVED = True +# params.output.spectra.HAS_TO_PLOT_SAVED = True +# params.output.spect_energy_budg.HAS_TO_PLOT_SAVED = True -params.output.periods_save.phys_fields = 1. +params.output.periods_save.phys_fields = 1.0 params.output.periods_save.spectra = 0.1 params.output.periods_save.spatial_means = 0.1 params.output.periods_save.spect_energy_budg = 0.5 @@ -51,16 +51,17 @@ params.preprocess.enable = True params.preprocess.forcing_const = 5e4 # params.preprocess.forcing_scale = 'unity' -params.preprocess.forcing_scale = 'enstrophy' +params.preprocess.forcing_scale = "enstrophy" params.preprocess.viscosity_const = 1.5 -params.preprocess.viscosity_scale = 'enstrophy_forcing' -params.preprocess.viscosity_type = 'hyper8' +params.preprocess.viscosity_scale = "enstrophy_forcing" +params.preprocess.viscosity_type = "hyper8" sim = Simul(params) sim.output.print_stdout( "Froude number ~ {:3e}, Viscosity = {:3e} + {:3e}".format( - (1./params.c2), sim.params.nu_2, sim.params.nu_8)) + (1.0 / params.c2), sim.params.nu_2, sim.params.nu_8 + ) +) sim.output.print_stdout(f"Forcing rate = {sim.params.forcing.forcing_rate}") sim.time_stepping.start() - diff --git a/scripts/launch/sw1l_no_flow.py b/scripts/launch/sw1l_no_flow.py index b035bc2c..f7c45687 100755 --- a/scripts/launch/sw1l_no_flow.py +++ b/scripts/launch/sw1l_no_flow.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -#coding=utf8 +# coding=utf8 # from math import pi @@ -8,26 +8,26 @@ params = Simul.create_default_params() -params.short_name_type_run = 'no_flow' -params.init_fields.type = 'constant' -params.init_fields.constant[...] = 0. -params.output.sub_directory = 'spect_energy_budg_tests_noflow' +params.short_name_type_run = "no_flow" +params.init_fields.type = "constant" +params.init_fields.constant[...] = 0.0 +params.output.sub_directory = "spect_energy_budg_tests_noflow" # --------Grid parameters--------- nh = 256 -Lh = 2. * pi +Lh = 2.0 * pi params.oper.nx = nh params.oper.ny = nh params.oper.Lx = Lh params.oper.Ly = Lh -delta_x = Lh/nh +delta_x = Lh / nh # -----Numerical method parameters---- -params.f = 0. -params.c2 = 10.**2 -params.oper.coef_dealiasing = 8./9 +params.f = 0.0 +params.c2 = 10.0 ** 2 +params.oper.coef_dealiasing = 8.0 / 9 params.time_stepping.USE_T_END = True -params.time_stepping.t_end = 5. +params.time_stepping.t_end = 5.0 params.time_stepping.it_end = 2000 params.time_stepping.USE_CFL = True params.time_stepping.deltat0 = 1e-3 @@ -35,19 +35,21 @@ # ----- Forcing parameters ----------- params.FORCING = True if params.FORCING: - params.forcing.type = 'waves' + params.forcing.type = "waves" params.forcing.nkmax_forcing = 8 params.forcing.nkmin_forcing = 5 - params.forcing.forcing_rate = 1. - k_max = pi / delta_x * params.oper.coef_dealiasing # Smallest resolved scale + params.forcing.forcing_rate = 1.0 + k_max = pi / delta_x * params.oper.coef_dealiasing # Smallest resolved scale k_d = k_max / 2.5 length_scale = pi / k_d - params.nu_8 = params.forcing.forcing_rate**(1./3) * length_scale ** (22./3) + params.nu_8 = params.forcing.forcing_rate ** (1.0 / 3) * length_scale ** ( + 22.0 / 3 + ) # params.nu_2 = params.forcing.forcing_rate ** (1./3) * length_scale ** (4./3) # ------Save file / StdOut parameters------------ params.output.periods_save.spatial_means = 0.1 -params.output.periods_save.phys_fields = 1. +params.output.periods_save.phys_fields = 1.0 params.output.periods_save.spectra = 0.1 params.output.periods_save.spect_energy_budg = 0.1 params.output.periods_save.increments = 0.5 @@ -64,17 +66,20 @@ # -------Preprocess parameters--------------- params.preprocess.enable = False -params.preprocess.init_field_scale = 'unity' -params.preprocess.forcing_const = 1. -params.preprocess.forcing_scale = 'unity' +params.preprocess.init_field_scale = "unity" +params.preprocess.forcing_const = 1.0 +params.preprocess.forcing_scale = "unity" params.preprocess.viscosity_const = 2.5 -params.preprocess.viscosity_scale = 'forcing' -params.preprocess.viscosity_type = 'hyper8' +params.preprocess.viscosity_scale = "forcing" +params.preprocess.viscosity_type = "hyper8" # -------Simulation commences------------- sim = Simul(params) if mpi.rank == 0: - print("Froude number ~ {:3e}, Viscosity = {:3e} + {:3e}".format( - (1./params.c2), sim.params.nu_2, sim.params.nu_8)) + print( + "Froude number ~ {:3e}, Viscosity = {:3e} + {:3e}".format( + (1.0 / params.c2), sim.params.nu_2, sim.params.nu_8 + ) + ) sim.time_stepping.start() diff --git a/scripts/ns2d.strat/coeff_diss.py b/scripts/ns2d.strat/coeff_diss.py index a0bc8e38..44c9b458 100644 --- a/scripts/ns2d.strat/coeff_diss.py +++ b/scripts/ns2d.strat/coeff_diss.py @@ -1,4 +1,3 @@ - """ coeff_diss.py ============= @@ -45,15 +44,18 @@ # SIMULATION gamma = args.gamma -F = np.sin(pi / 4) # F = omega_l / N -sigma = 1 # sigma = omega_l / (pi * f_cf); f_cf freq time correlation forcing in s-1 +F = np.sin(pi / 4) # F = omega_l / N +sigma = ( + 1 # sigma = omega_l / (pi * f_cf); f_cf freq time correlation forcing in s-1 +) nu_8 = 1e-16 NO_SHEAR_MODES = False -t_end = 8. +t_end = 8.0 PLOT_FIGURES = False ################### + def load_mean_spect_energy_budg(sim, tmin=0, tmax=1000): """ Loads data spect_energy_budget. @@ -73,41 +75,50 @@ def load_mean_spect_energy_budg(sim, tmin=0, tmax=1000): imax_plot = np.argmin(abs(times - tmax)) dset_dissE_kx = dset_dissEK_kx + dset_dissEA_kx - dissE_kx = dset_dissE_kx[imin_plot:imax_plot + 1].mean(0) + dissE_kx = dset_dissE_kx[imin_plot : imax_plot + 1].mean(0) dset_dissE_ky = dset_dissEK_ky + dset_dissEA_ky - dissE_ky = dset_dissE_ky[imin_plot:imax_plot + 1].mean(0) + dissE_ky = dset_dissE_ky[imin_plot : imax_plot + 1].mean(0) return kxE, kyE, dissE_kx, dissE_ky + def get_state_from_sim(sim): """Returns the state from a simulation.""" # Take the state - b_fft = sim.state.get_var('b_fft') - rot_fft = sim.state.get_var('rot_fft') + b_fft = sim.state.get_var("b_fft") + rot_fft = sim.state.get_var("rot_fft") return rot_fft, b_fft + def compute_diff(idx_dealiasing, idy_dealiasing, dissE_kx, dissE_ky): idx_diss_max = np.argmax(abs(dissE_kx)) idy_diss_max = np.argmax(abs(dissE_ky)) # Computes difference - diff_x = idx_dealiasing - idx_diss_max - diff_y = idy_dealiasing - idy_diss_max + diff_x = idx_dealiasing - idx_diss_max + diff_y = idy_dealiasing - idy_diss_max diff = np.argmax([diff_x, diff_y]) if diff == 0: - print("diff_x = idx_dealiasing - idx_diss_max", idx_dealiasing - idx_diss_max) - diff = idx_dealiasing - idx_diss_max + print( + "diff_x = idx_dealiasing - idx_diss_max", + idx_dealiasing - idx_diss_max, + ) + diff = idx_dealiasing - idx_diss_max else: - print("diff_y = idy_dealiasing - idy_diss_max", idy_dealiasing - idy_diss_max) - diff = idy_dealiasing - idy_diss_max + print( + "diff_y = idy_dealiasing - idy_diss_max", + idy_dealiasing - idy_diss_max, + ) + diff = idy_dealiasing - idy_diss_max return diff, idx_diss_max, idy_diss_max + def normalization_initialized_field(sim, coef_norm=1e-4): """Normalizes the initialized field. (ONLY if nx != ny)""" if sim.params.oper.nx != sim.params.oper.ny: @@ -116,19 +127,19 @@ def normalization_initialized_field(sim, coef_norm=1e-4): raise ValueError("sim.params.forcing.key_forced should be ap_fft.") KX = sim.oper.KX - cond = KX == 0. + cond = KX == 0.0 - ux_fft = sim.state.get_var('ux_fft') - uy_fft = sim.state.get_var('uy_fft') - b_fft = sim.state.get_var('b_fft') + ux_fft = sim.state.get_var("ux_fft") + uy_fft = sim.state.get_var("uy_fft") + b_fft = sim.state.get_var("b_fft") - ux_fft[cond] = 0. - uy_fft[cond] = 0. - b_fft[cond] = 0. + ux_fft[cond] = 0.0 + uy_fft[cond] = 0.0 + b_fft[cond] = 0.0 # Compute energy after ux_fft[kx=0] uy_fft[kx=0] b_fft[kx=0] - ek_fft = (np.abs(ux_fft)**2 + np.abs(uy_fft)**2)/2 - ea_fft = ((np.abs(b_fft)/params.N)**2)/2 + ek_fft = (np.abs(ux_fft) ** 2 + np.abs(uy_fft) ** 2) / 2 + ea_fft = ((np.abs(b_fft) / params.N) ** 2) / 2 e_fft = ek_fft + ea_fft energy_before_norm = sim.output.sum_wavenumbers(e_fft) @@ -139,8 +150,10 @@ def normalization_initialized_field(sim, coef_norm=1e-4): nkmax_forcing = params.forcing.nkmax_forcing nkmin_forcing = params.forcing.nkmin_forcing - k_f = ((nkmax_forcing + nkmin_forcing) / 2) * max(2 * pi / Lx, 2 * pi / Lz) - energy_f = params.forcing.forcing_rate**(2/7) * (2 * pi / k_f)**7 + k_f = ((nkmax_forcing + nkmin_forcing) / 2) * max( + 2 * pi / Lx, 2 * pi / Lz + ) + energy_f = params.forcing.forcing_rate ** (2 / 7) * (2 * pi / k_f) ** 7 coef = np.sqrt(coef_norm * energy_f / energy_before_norm) @@ -156,18 +169,21 @@ def normalization_initialized_field(sim, coef_norm=1e-4): pass return sim + def _compute_ikdiss(dissE_kx, dissE_ky): idx_diss_max = np.argmax(abs(dissE_kx)) idy_diss_max = np.argmax(abs(dissE_ky)) return idx_diss_max, idy_diss_max + def _compute_ikmax(kxE, kyE): idx_dealiasing = np.argmin(abs(kxE - sim.oper.kxmax_dealiasing)) idy_dealiasing = np.argmin(abs(kyE - sim.oper.kymax_dealiasing)) return idx_dealiasing, idy_dealiasing + def compute_delta_ik(kxE, kyE, dissE_kx, dissE_ky): idx_diss_max, idy_diss_max = _compute_ikdiss(dissE_kx, dissE_ky) idx_dealiasing, idy_dealiasing = _compute_ikmax(kxE, kyE) @@ -180,6 +196,7 @@ def compute_delta_ik(kxE, kyE, dissE_kx, dissE_ky): return idx_target, idy_target, delta_ikx, delta_iky + def compute_energy_spatial(sim): """ Compute energy without energy in shear modes """ dict_spatial = sim.output.spatial_means.load() @@ -191,6 +208,7 @@ def compute_energy_spatial(sim): P_tot = PK_tot + PA_tot return E, t, P_tot + def modify_factor(sim): params = _deepcopy(sim.params) @@ -199,7 +217,7 @@ def modify_factor(sim): sim_old = sim params.nu_8 = params.nu_8 * factor - params.init_fields.type = 'in_script' + params.init_fields.type = "in_script" params.time_stepping.t_end = t_end return params, nu_8_old @@ -209,6 +227,7 @@ def write_to_file(path, to_print, mode="a"): with open(path, mode) as f: f.write(to_print) + def make_float_value_for_path(value): value_not = str(value) if "." in value_not: @@ -216,16 +235,17 @@ def make_float_value_for_path(value): else: return value_not + def check_dissipation(): ratio_x = dissE_kx[idx_diss_max] / dissE_kx[idx_dealiasing - 1] ratio_y = dissE_ky[idy_diss_max] / dissE_ky[idy_dealiasing - 1] - idx_target, idy_target, delta_ikx, delta_iky = compute_delta_ik(kxE, kyE, dissE_kx, dissE_ky) + idx_target, idy_target, delta_ikx, delta_iky = compute_delta_ik( + kxE, kyE, dissE_kx, dissE_ky + ) if time_total > 1000: - print( - "The stationarity has not " + \ - f"reached after {it} simulations.") + print("The stationarity has not " + f"reached after {it} simulations.") should_I_stop = "non_stationarity" if ratio_x > threshold_ratio and ratio_y > threshold_ratio: @@ -248,14 +268,15 @@ def check_dissipation(): print(f"Checking stationarity... with nu8 = {nu_8_old}") E, t, P_tot = compute_energy_spatial(sim) ratio = np.mean(np.diff(E[2:]) / np.diff(t[2:])) - if (ratio / injection_energy_0) < 0.5 and \ - abs(nu_8_old - params.nu_8) / params.nu_8 < 0.05: + if (ratio / injection_energy_0) < 0.5 and abs( + nu_8_old - params.nu_8 + ) / params.nu_8 < 0.05: print(f"Stationarity is reached.\n nu_8 = {params.nu_8}") - factor = 1. + factor = 1.0 should_I_stop = True else: should_I_stop = False - factor = 1. + factor = 1.0 else: factor = 1 + (1 / min(ratio_x, ratio_y)) @@ -263,6 +284,7 @@ def check_dissipation(): return factor, should_I_stop + ################################################# gamma_not = make_float_value_for_path(gamma) @@ -291,7 +313,9 @@ def check_dissipation(): if len(paths_sim) == 0: # Make FIRST SIMULATION - params = make_parameters_simulation(gamma, F, sigma, nu_8, t_end=t_end, NO_SHEAR_MODES=NO_SHEAR_MODES) + params = make_parameters_simulation( + gamma, F, sigma, nu_8, t_end=t_end, NO_SHEAR_MODES=NO_SHEAR_MODES + ) sim = Simul(params) sim = normalization_initialized_field(sim) sim.time_stepping.start() @@ -331,7 +355,7 @@ def check_dissipation(): time_total += sim.time_stepping.t energy, t, P_tot = compute_energy_spatial(sim) -energy = np.mean(energy[len(energy) // 2:]) +energy = np.mean(energy[len(energy) // 2 :]) injection_energy_0 = P_tot[2] energies = [] @@ -341,20 +365,25 @@ def check_dissipation(): # Write results into temporal file if mpi.rank == 0: - to_print = ("####\n" - "t = {:.4e} \n" - "E = {:.4e} \n" - "nu8 = {:.4e} \n" - "factor = {:.4e} \n").format( - time_total, energy, sim.params.nu_8, 1) + to_print = ( + "####\n" + "t = {:.4e} \n" + "E = {:.4e} \n" + "nu8 = {:.4e} \n" + "factor = {:.4e} \n" + ).format(time_total, energy, sim.params.nu_8, 1) write_to_file(path_file_write2, to_print, mode="w") # Compute the data spectra energy budget -kxE, kyE, dissE_kx, dissE_ky = load_mean_spect_energy_budg(sim, tmin=2., tmax=1000) +kxE, kyE, dissE_kx, dissE_ky = load_mean_spect_energy_budg( + sim, tmin=2.0, tmax=1000 +) idx_diss_max, idy_diss_max = _compute_ikdiss(dissE_kx, dissE_ky) -idx_dealiasing, idy_dealiasing = _compute_ikmax(kxE, kyE) -idx_target, idy_target, delta_ikx, delta_iky = compute_delta_ik(kxE, kyE, dissE_kx, dissE_ky) +idx_dealiasing, idy_dealiasing = _compute_ikmax(kxE, kyE) +idx_target, idy_target, delta_ikx, delta_iky = compute_delta_ik( + kxE, kyE, dissE_kx, dissE_ky +) # Plot dissipation if PLOT_FIGURES: @@ -363,18 +392,28 @@ def check_dissipation(): ax.set_title("$D(k_y)$") ax.set_xlabel("$k_y$") ax.set_ylabel("$D(k_y)$") - ax.plot(kyE, dissE_ky, label="nu8 = {:.2e}, diff = {}".format( - sim.params.nu_8, abs(idy_diss_max - idy_dealiasing))) - ax.plot(sim.oper.kymax_dealiasing, 0, 'xr') + ax.plot( + kyE, + dissE_ky, + label="nu8 = {:.2e}, diff = {}".format( + sim.params.nu_8, abs(idy_diss_max - idy_dealiasing) + ), + ) + ax.plot(sim.oper.kymax_dealiasing, 0, "xr") ax.axvline(x=sim.oper.deltaky * idy_target, color="k") fig2, ax2 = plt.subplots() ax2.set_title("$D(k_x)$") ax2.set_xlabel("$k_x$") ax2.set_ylabel("$D(k_x)$") - ax2.plot(kxE, dissE_kx, label="nu8 = {:.2e}, diff = {}".format( - sim.params.nu_8, abs(idx_diss_max - idx_dealiasing))) - ax2.plot(sim.oper.kxmax_dealiasing, 0, 'xr') + ax2.plot( + kxE, + dissE_kx, + label="nu8 = {:.2e}, diff = {}".format( + sim.params.nu_8, abs(idx_diss_max - idx_dealiasing) + ), + ) + ax2.plot(sim.oper.kxmax_dealiasing, 0, "xr") ax2.axvline(x=sim.oper.deltakx * idx_target, color="k") ax.legend() @@ -390,17 +429,17 @@ def check_dissipation(): ax3.set_xlabel("times") ax3.set_ylabel(r"$\nu_8$") - ax3.plot(time_total, viscosities[-1], '.') + ax3.plot(time_total, viscosities[-1], ".") # Energy Vs time fig4, ax4 = plt.subplots() - ax4.plot(time_total, energy, '.') + ax4.plot(time_total, energy, ".") ax4.set_xlabel("times") ax4.set_ylabel("Energy") # Factor Vs time fig5, ax5 = plt.subplots() - ax5.plot(time_total, 1, '.') + ax5.plot(time_total, 1, ".") ax5.set_xlabel("times") ax5.set_ylabel("Factor") @@ -444,31 +483,43 @@ def check_dissipation(): # Add values to time array and energy array time_total += sim.time_stepping.t energy, t, P_tot = compute_energy_spatial(sim) - energy = np.mean(energy[len(energy) // 2:]) + energy = np.mean(energy[len(energy) // 2 :]) energies.append(energy) viscosities.append(sim.params.nu_8) # Computes new index k_max_dissipation kxE, kyE, dissE_kx, dissE_ky = load_mean_spect_energy_budg( - sim, tmin=2., tmax=1000) + sim, tmin=2.0, tmax=1000 + ) idx_diss_max, idy_diss_max = _compute_ikdiss(dissE_kx, dissE_ky) - idx_dealiasing, idy_dealiasing = _compute_ikmax(kxE, kyE) + idx_dealiasing, idy_dealiasing = _compute_ikmax(kxE, kyE) # Write results into temporal file - to_print = ("####\n" - "t = {:.4e} \n" - "E = {:.4e} \n" - "nu8 = {:.4e} \n" - "factor = {:.4e} \n").format( - time_total, energy, sim.params.nu_8, factor) + to_print = ( + "####\n" + "t = {:.4e} \n" + "E = {:.4e} \n" + "nu8 = {:.4e} \n" + "factor = {:.4e} \n" + ).format(time_total, energy, sim.params.nu_8, factor) write_to_file(path_file_write2, to_print, mode="a") if PLOT_FIGURES: - ax.plot(kyE, dissE_ky, label="nu8 = {:.2e}, diff = {}".format( - params.nu_8, abs(idy_diss_max - idy_dealiasing))) - ax2.plot(kxE, dissE_kx, label="nu8 = {:.2e}, diff = {}".format( - params.nu_8, abs(idx_diss_max - idx_dealiasing))) + ax.plot( + kyE, + dissE_ky, + label="nu8 = {:.2e}, diff = {}".format( + params.nu_8, abs(idy_diss_max - idy_dealiasing) + ), + ) + ax2.plot( + kxE, + dissE_kx, + label="nu8 = {:.2e}, diff = {}".format( + params.nu_8, abs(idx_diss_max - idx_dealiasing) + ), + ) fig.canvas.draw() fig2.canvas.draw() @@ -483,7 +534,7 @@ def check_dissipation(): ax4.autoscale() fig4.canvas.draw() - ax5.plot(time_total, factor, 'x') + ax5.plot(time_total, factor, "x") ax5.autoscale() fig5.canvas.draw() @@ -491,14 +542,12 @@ def check_dissipation(): if mpi.rank == 0: if len(paths_sim) == 0: - to_print = ("gamma,nx,nu8 \n") - to_print += ("{},{},{} \n".format( - gamma, sim.params.oper.nx, params.nu_8)) + to_print = "gamma,nx,nu8 \n" + to_print += "{},{},{} \n".format(gamma, sim.params.oper.nx, params.nu_8) mode_write = "w" else: - to_print = ("{},{},{} \n".format( - gamma, sim.params.oper.nx, params.nu_8)) + to_print = "{},{},{} \n".format(gamma, sim.params.oper.nx, params.nu_8) mode_write = "a" write_to_file(path_file_write, to_print, mode=mode_write) diff --git a/scripts/ns2d.strat/compute_anisotropy.py b/scripts/ns2d.strat/compute_anisotropy.py index 20009b0b..384b4f4d 100644 --- a/scripts/ns2d.strat/compute_anisotropy.py +++ b/scripts/ns2d.strat/compute_anisotropy.py @@ -13,6 +13,7 @@ from fluidsim import load_params_simul + def _compute_array_times_from_path(path_simulation): """ Compute array with times from path simulation. @@ -28,7 +29,8 @@ def _compute_array_times_from_path(path_simulation): for path in paths_phys_files: if not "=" in path.split("state_phys_t")[1].split(".nc")[0]: times_phys_files.append( - float(path.split("state_phys_t")[1].split(".nc")[0])) + float(path.split("state_phys_t")[1].split(".nc")[0]) + ) else: continue return np.asarray(times_phys_files) @@ -83,14 +85,17 @@ def compute_anisotropy(path_simulation, tmin=None): delta_kx = 2 * np.pi / params.oper.Lx # Compute spatial averaged energy from spectra - EK_ux = np.sum(np.mean(spectrumkx_EK_ux[itmin:,:], axis=0) * delta_kx) - EK = np.sum(np.mean(spectrumkx_EK[itmin:,:], axis=0) * delta_kx) + EK_ux = np.sum(np.mean(spectrumkx_EK_ux[itmin:, :], axis=0) * delta_kx) + EK = np.sum(np.mean(spectrumkx_EK[itmin:, :], axis=0) * delta_kx) return EK_ux / EK + if __name__ == "__main__": - path_simulation = ("/fsnet/project/meige/2015/15DELDUCA/DataSim/" + - "sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma05_2018-08-14_10-01-22") + path_simulation = ( + "/fsnet/project/meige/2015/15DELDUCA/DataSim/" + + "sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma05_2018-08-14_10-01-22" + ) anisotropy = compute_anisotropy(path_simulation) print(f"anisotropy = {anisotropy}") diff --git a/scripts/ns2d.strat/compute_flow_features.py b/scripts/ns2d.strat/compute_flow_features.py index 40827d1d..0b62fdeb 100644 --- a/scripts/ns2d.strat/compute_flow_features.py +++ b/scripts/ns2d.strat/compute_flow_features.py @@ -15,7 +15,7 @@ # Argparse arguments nx = 3840 -n_files_average = 50 # Number of files to perform the time average +n_files_average = 50 # Number of files to perform the time average # Create paths path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" @@ -50,31 +50,36 @@ with h5py.File(path_file, "r") as f: ux = f["state_phys"]["ux"][...] uz = f["state_phys"]["uy"][...] - anisotropies.append(np.mean(ux**2) / np.mean(uz**2)) + anisotropies.append(np.mean(ux ** 2) / np.mean(uz ** 2)) anisotropies_gammas.append(np.mean(anisotropies)) # Compute ratio D(k_x)/epsilon print("Computing ratio dissipation for gamma {}...".format(gammas[ipath])) with h5py.File(path + "/spect_energy_budg.h5", "r") as f: - kx = f['kxE'][...] - kz = f['kyE'][...] - dset_dissEKu_kx = f['dissEKu_kx'] - dset_dissEKv_kx = f['dissEKv_kx'] - dset_dissEA_kx = f['dissEA_kx'] + kx = f["kxE"][...] + kz = f["kyE"][...] + dset_dissEKu_kx = f["dissEKu_kx"] + dset_dissEKv_kx = f["dissEKv_kx"] + dset_dissEA_kx = f["dissEA_kx"] delta_kx = kx[1] - kx[0] delta_kz = kz[1] - kz[0] - dissEK_kx = (dset_dissEKu_kx[-n_files_average:] + \ - dset_dissEKv_kx[-n_files_average:]) + dissEK_kx = ( + dset_dissEKu_kx[-n_files_average:] + + dset_dissEKv_kx[-n_files_average:] + ) dissEA_kx = dset_dissEA_kx[-n_files_average:] dissE_kx = (dissEK_kx + dissEA_kx).mean(0) D_kx = dissE_kx.cumsum() * delta_kx # Compute k_fx - k_fx = (np.sin(params.forcing.tcrandom_anisotropic.angle) * - params.forcing.nkmax_forcing * max(delta_kx, delta_kz)) - ik_fx = np.argmin(abs(kx - k_fx )) + k_fx = ( + np.sin(params.forcing.tcrandom_anisotropic.angle) + * params.forcing.nkmax_forcing + * max(delta_kx, delta_kz) + ) + ik_fx = np.argmin(abs(kx - k_fx)) # Compute ratio ratio_dissipations.append(D_kx[ik_fx] / D_kx[-1]) @@ -82,11 +87,11 @@ fig1, ax1 = plt.subplots() ax1.set_xlabel(r"$\gamma$") ax1.set_ylabel(r"$U_x^2/U_z^2$") -ax1.plot(gammas, anisotropies_gammas, 'ro') +ax1.plot(gammas, anisotropies_gammas, "ro") fig2, ax2 = plt.subplots() ax2.set_xlabel(r"$\gamma$") ax2.set_ylabel(r"$D(k_{fx})/D(k_{x, max})$") -ax2.plot(gammas, ratio_dissipations, 'bo') +ax2.plot(gammas, ratio_dissipations, "bo") plt.show() diff --git a/scripts/ns2d.strat/compute_length_scales.py b/scripts/ns2d.strat/compute_length_scales.py index 566279fa..2499a7c4 100644 --- a/scripts/ns2d.strat/compute_length_scales.py +++ b/scripts/ns2d.strat/compute_length_scales.py @@ -10,6 +10,7 @@ from fluidsim import load_params_simul + def compute_length_scales(path_simulation, tmin=None): """ Compute length scales. @@ -61,16 +62,17 @@ def compute_length_scales(path_simulation, tmin=None): delta_ky = np.median(np.diff(abs(ky))) # Compute vertical length scale - lx = (np.sum(spectrum1Dkx_EK_ux * delta_kx) / - np.sum(kx * spectrum1Dkx_EK_ux * delta_kx)) + lx = np.sum(spectrum1Dkx_EK_ux * delta_kx) / np.sum( + kx * spectrum1Dkx_EK_ux * delta_kx + ) - lz = (np.sum(spectrum1Dky_EK_uy * delta_ky) / - np.sum(ky * spectrum1Dky_EK_uy * delta_ky)) + lz = np.sum(spectrum1Dky_EK_uy * delta_ky) / np.sum( + ky * spectrum1Dky_EK_uy * delta_ky + ) return lx, lz - if __name__ == "__main__": # path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma02_2018-08-14_09-59-55" diff --git a/scripts/ns2d.strat/compute_ratio_dissipation.py b/scripts/ns2d.strat/compute_ratio_dissipation.py index 303208bb..a4b9a1ea 100644 --- a/scripts/ns2d.strat/compute_ratio_dissipation.py +++ b/scripts/ns2d.strat/compute_ratio_dissipation.py @@ -9,6 +9,7 @@ from fluidsim import load_params_simul + def compute_ratio_dissipation(path_simulation, tmin=None): """ Compute ratio dissipation from path simulation. @@ -29,9 +30,9 @@ def compute_ratio_dissipation(path_simulation, tmin=None): with h5py.File(path_simulation + "/spect_energy_budg.h5", "r") as f: times = f["times"][...] kx = f["kxE"][...] - kz = f['kyE'][...] - dset_dissEKu_kx = f['dissEKu_kx'] - dset_dissEKv_kx = f['dissEKv_kx'] + kz = f["kyE"][...] + dset_dissEKu_kx = f["dissEKu_kx"] + dset_dissEKv_kx = f["dissEKv_kx"] dset_dissEA_kx = f["dissEA_kx"] # Compute itmin time average @@ -50,8 +51,11 @@ def compute_ratio_dissipation(path_simulation, tmin=None): D_kx = dissE_kx.cumsum() * delta_kx # Compute k_fx - k_fx = (np.sin(params.forcing.tcrandom_anisotropic.angle) * - params.forcing.nkmax_forcing * max(delta_kx, delta_kz)) + k_fx = ( + np.sin(params.forcing.tcrandom_anisotropic.angle) + * params.forcing.nkmax_forcing + * max(delta_kx, delta_kz) + ) ik_fx = np.argmin(abs(kx - k_fx)) ratio_dissipation_old = D_kx[ik_fx] / D_kx[-1] @@ -62,10 +66,13 @@ def compute_ratio_dissipation(path_simulation, tmin=None): # return D_kx[ik_half] / D_kx[ik_fx] return kx[ik_half] / k_fx + if __name__ == "__main__": - path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/" + \ - "sim960_no_shear_modes/" + \ - "NS2D.strat_960x240_S2pix1.571_F07_gamma02_2018-08-10_14-45-06" + path_simulation = ( + "/fsnet/project/meige/2015/15DELDUCA/DataSim/" + + "sim960_no_shear_modes/" + + "NS2D.strat_960x240_S2pix1.571_F07_gamma02_2018-08-10_14-45-06" + ) ratio = compute_ratio_dissipation(path_simulation, tmin=None) print(f"Ratio between ik_half / ik_fx = {ratio}") diff --git a/scripts/ns2d.strat/compute_reynolds_froude.py b/scripts/ns2d.strat/compute_reynolds_froude.py index 8d0e6f80..3681802d 100644 --- a/scripts/ns2d.strat/compute_reynolds_froude.py +++ b/scripts/ns2d.strat/compute_reynolds_froude.py @@ -10,6 +10,7 @@ from fluidsim import load_params_simul + def _compute_epsilon_from_path(path_simulation, tmin=None): """ Computes the mean dissipation from tmin @@ -32,13 +33,13 @@ def _compute_epsilon_from_path(path_simulation, tmin=None): epsK = np.empty(nt) for il in range(nt): - line = lines_t[il] - words = line.split() - t[il] = float(words[2]) + line = lines_t[il] + words = line.split() + t[il] = float(words[2]) - line = lines_epsK[il] - words = line.split() - epsK[il] = float(words[2]) + line = lines_epsK[il] + words = line.split() + epsK[il] = float(words[2]) # Compute start time average itmin dt_spatial = np.median(np.diff(t)) @@ -51,6 +52,7 @@ def _compute_epsilon_from_path(path_simulation, tmin=None): return np.mean(epsK[itmin:], axis=0) + def _compute_lx_from_path(path_simulation, tmin=None): """ Compute horizontal length from path using appendix B. Brethouwer 2007 @@ -81,8 +83,10 @@ def _compute_lx_from_path(path_simulation, tmin=None): delta_kx = np.median(np.diff(kx)) # Compute horizontal length scale Brethouwer 2007 - return (np.sum(spectrum1Dkx_EK_ux * delta_kx) / - np.sum(kx * spectrum1Dkx_EK_ux * delta_kx)) + return np.sum(spectrum1Dkx_EK_ux * delta_kx) / np.sum( + kx * spectrum1Dkx_EK_ux * delta_kx + ) + def compute_buoyancy_reynolds(path_simulation, tmin=None): """ @@ -92,14 +96,15 @@ def compute_buoyancy_reynolds(path_simulation, tmin=None): epsK = _compute_epsilon_from_path(path_simulation, tmin=tmin) lx = _compute_lx_from_path(path_simulation, tmin=tmin) - F_h = ((epsK / lx**2)**(1/3)) * (1 / params.N) + F_h = ((epsK / lx ** 2) ** (1 / 3)) * (1 / params.N) - eta_8 = (params.nu_8 ** 3 / epsK)**(1/22) - Re_8 = (lx/eta_8)**(22/3) + eta_8 = (params.nu_8 ** 3 / epsK) ** (1 / 22) + Re_8 = (lx / eta_8) ** (22 / 3) - R_b = Re_8 * F_h**8 + R_b = Re_8 * F_h ** 8 return F_h, Re_8, R_b + # path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma02_2018-08-14_09-59-55" # F_h, Re_8, R_b = compute_buoyancy_reynolds(path_simulation) # print("F_h", F_h) diff --git a/scripts/ns2d.strat/compute_spectrum_kykx.py b/scripts/ns2d.strat/compute_spectrum_kykx.py index ebf37b8b..41ea21b8 100644 --- a/scripts/ns2d.strat/compute_spectrum_kykx.py +++ b/scripts/ns2d.strat/compute_spectrum_kykx.py @@ -71,7 +71,7 @@ kz = 2 * pi * np.fft.fftfreq(nz, Lz / nz) KX, KZ = np.meshgrid(kx, kz) -omega_k = sim.params.N * (KX / np.sqrt(KX**2 + KZ**2)) +omega_k = sim.params.N * (KX / np.sqrt(KX ** 2 + KZ ** 2)) # Create 3D_arrays ux_fft_arr = np.empty([itmax - itmin, nz, nx], dtype="complex") @@ -90,8 +90,12 @@ ux_fft_arr[ifile, :, :] = np.fft.fft2(ux) uz_fft_arr[ifile, :, :] = np.fft.fft2(uz) b_fft_arr[ifile, :, :] = np.fft.fft2(b) - ap_fft_arr[ifile, :, :] = N**2 * np.fft.fft2(uz) + 1j * omega_k * np.fft.fft2(b) - am_fft_arr[ifile, :, :] = N**2 * np.fft.fft2(uz) - 1j * omega_k * np.fft.fft2(b) + ap_fft_arr[ifile, :, :] = N ** 2 * np.fft.fft2( + uz + ) + 1j * omega_k * np.fft.fft2(b) + am_fft_arr[ifile, :, :] = N ** 2 * np.fft.fft2( + uz + ) - 1j * omega_k * np.fft.fft2(b) # Time average ux_fft_arr = np.mean(ux_fft_arr, axis=0) @@ -122,16 +126,15 @@ # data = abs(am_fft_arr)**2 if key == "ap": - data = np.log10(abs(ap_fft_arr)**2) + data = np.log10(abs(ap_fft_arr) ** 2) text = r"$\hat{a}_+$" elif key == "am": - data = np.log10(abs(am_fft_arr)**2) + data = np.log10(abs(am_fft_arr) ** 2) text = r"$\hat{a}_-$" else: raise ValueError("Not implemented") - ### Data ikx = np.argmin(abs(kx - 200)) ikz = np.argmin(abs(kz - 148)) @@ -151,30 +154,21 @@ # ) ax1.pcolormesh( - KX[0:ikz, 0:ikx], - KZ[0:ikz, 0:ikx], - data[0:ikz, 0:ikx], - vmin=5, - vmax=8) + KX[0:ikz, 0:ikx], KZ[0:ikz, 0:ikx], data[0:ikz, 0:ikx], vmin=5, vmax=8 +) # # Plot kx > 0 and kz < 0 -KX_grid2 = KX[ikz_negative - 1:, 0:ikx] +KX_grid2 = KX[ikz_negative - 1 :, 0:ikx] KZ_grid2 = np.empty_like(KX_grid2) -KZ_grid2[0:-1,:] = KZ[ikz_negative:, 0:ikx] -KZ_grid2[-1,:] = KZ[0, 0:ikx] +KZ_grid2[0:-1, :] = KZ[ikz_negative:, 0:ikx] +KZ_grid2[-1, :] = KZ[0, 0:ikx] data_grid2 = np.empty_like(KX_grid2) -data_grid2[0:-1,:] = data[ikz_negative:, 0:ikx] -data_grid2[-1,:] = data[0, 0:ikx] +data_grid2[0:-1, :] = data[ikz_negative:, 0:ikx] +data_grid2[-1, :] = data[0, 0:ikx] -ax1.pcolormesh( - KX_grid2, - KZ_grid2, - data_grid2, - vmin=5, - vmax=8 -) +ax1.pcolormesh(KX_grid2, KZ_grid2, data_grid2, vmin=5, vmax=8) ############## @@ -213,24 +207,47 @@ # Create a Rectangle patch deltak = max(sim.oper.deltakx, sim.oper.deltaky) -x_rect = np.sin( - sim.params.forcing.tcrandom_anisotropic.angle) * deltak * sim.params.forcing.nkmin_forcing +x_rect = ( + np.sin(sim.params.forcing.tcrandom_anisotropic.angle) + * deltak + * sim.params.forcing.nkmin_forcing +) -z_rect = np.cos( - sim.params.forcing.tcrandom_anisotropic.angle) * deltak * sim.params.forcing.nkmin_forcing +z_rect = ( + np.cos(sim.params.forcing.tcrandom_anisotropic.angle) + * deltak + * sim.params.forcing.nkmin_forcing +) width = abs( - x_rect - np.sin(sim.params.forcing.tcrandom_anisotropic.angle) * deltak * sim.params.forcing.nkmax_forcing) + x_rect + - np.sin(sim.params.forcing.tcrandom_anisotropic.angle) + * deltak + * sim.params.forcing.nkmax_forcing +) height = abs( - z_rect - np.cos(sim.params.forcing.tcrandom_anisotropic.angle) * deltak * sim.params.forcing.nkmax_forcing) + z_rect + - np.cos(sim.params.forcing.tcrandom_anisotropic.angle) + * deltak + * sim.params.forcing.nkmax_forcing +) -rect1 = patches.Rectangle((x_rect,z_rect),width,height,linewidth=1,edgecolor='r',facecolor='none') +rect1 = patches.Rectangle( + (x_rect, z_rect), width, height, linewidth=1, edgecolor="r", facecolor="none" +) ax1.add_patch(rect1) if sim.params.forcing.tcrandom_anisotropic.kz_negative_enable: - rect2 = patches.Rectangle((x_rect,-(z_rect + height)),width,height,linewidth=1,edgecolor='r',facecolor='none') + rect2 = patches.Rectangle( + (x_rect, -(z_rect + height)), + width, + height, + linewidth=1, + edgecolor="r", + facecolor="none", + ) ax1.add_patch(rect2) @@ -241,10 +258,10 @@ width=2 * sim.params.forcing.nkmin_forcing * deltak, height=2 * sim.params.forcing.nkmin_forcing * deltak, angle=0, - theta1=-90., - theta2=90., + theta1=-90.0, + theta2=90.0, linestyle="-.", - color="red" + color="red", ) ) ax1.add_patch( @@ -254,9 +271,9 @@ height=2 * sim.params.forcing.nkmax_forcing * deltak, angle=0, theta1=-90, - theta2=90., + theta2=90.0, linestyle="-.", - color="red" + color="red", ) ) diff --git a/scripts/ns2d.strat/compute_vertical_length.py b/scripts/ns2d.strat/compute_vertical_length.py index 27e66509..2017df29 100644 --- a/scripts/ns2d.strat/compute_vertical_length.py +++ b/scripts/ns2d.strat/compute_vertical_length.py @@ -10,6 +10,7 @@ from fluidsim import load_params_simul + def compute_length_scales(path_simulation, tmin=None): """ Compute length scales. @@ -63,16 +64,17 @@ def compute_length_scales(path_simulation, tmin=None): delta_ky = np.median(np.diff(abs(ky))) # Compute vertical length scale - lx = (np.sum(spectrum1Dkx_EK_ux * delta_kx) / - np.sum(kx * spectrum1Dkx_EK_ux * delta_kx)) + lx = np.sum(spectrum1Dkx_EK_ux * delta_kx) / np.sum( + kx * spectrum1Dkx_EK_ux * delta_kx + ) - lz = (np.sum(spectrum1Dky_EK_uy * delta_ky) / - np.sum(ky * spectrum1Dky_EK_uy * delta_ky)) + lz = np.sum(spectrum1Dky_EK_uy * delta_ky) / np.sum( + ky * spectrum1Dky_EK_uy * delta_ky + ) return lx, lz - if __name__ == "__main__": # path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma02_2018-08-14_09-59-55" diff --git a/scripts/ns2d.strat/flow_features.py b/scripts/ns2d.strat/flow_features.py index e09b8c85..8dce93cd 100644 --- a/scripts/ns2d.strat/flow_features.py +++ b/scripts/ns2d.strat/flow_features.py @@ -32,7 +32,8 @@ def compute_flow_features_from_sim(path_simulation, SAVE=False): """ if SAVE and os.path.exists( - os.path.join(path_simulation, "flow_features.txt")): + os.path.join(path_simulation, "flow_features.txt") + ): pass else: # Computes all the flow features @@ -43,13 +44,15 @@ def compute_flow_features_from_sim(path_simulation, SAVE=False): if SAVE: path_save = os.path.join(path_simulation, "flow_features.txt") - to_print = (f"anisotropy = {anisotropy} \n" + - f"ratio_diss = {ratio_diss} \n" + - f"F_h = {F_h} \n" + - f"Re_8 = {Re_8} \n" + - f"R_b = {R_b} \n" + - f"l_x = {l_x} \n" + - f"l_z = {l_z} \n") + to_print = ( + f"anisotropy = {anisotropy} \n" + + f"ratio_diss = {ratio_diss} \n" + + f"F_h = {F_h} \n" + + f"Re_8 = {Re_8} \n" + + f"R_b = {R_b} \n" + + f"l_x = {l_x} \n" + + f"l_z = {l_z} \n" + ) # Checks if the file flow_features.txt exists if os.path.exists(path_save): @@ -59,6 +62,7 @@ def compute_flow_features_from_sim(path_simulation, SAVE=False): f.write(to_print) return anisotropy, ratio_diss, F_h, Re_8, R_b, l_x, l_z + def get_features_from_sim(path_simulation): """ Gets features of a simulation @@ -87,12 +91,13 @@ def get_features_from_sim(path_simulation): if line.startswith("l_z = "): l_z = float(line.split()[2]) - return anisotropy, ratio_diss, F_h, Re_8, R_b, l_x, l_z + def _get_resolution_from_dir(path_simulation): return path_simulation.split("NS2D.strat_")[1].split("x")[0] + if __name__ == "__main__": # path_simulation = ("/fsnet/project/meige/2015/15DELDUCA/DataSim/" + # "sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma05_2018-08-14_10-01-22") @@ -103,18 +108,21 @@ def _get_resolution_from_dir(path_simulation): # Compute flow features from all simulations... path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" - directories = ["sim960_no_shear_modes", - "sim960_no_shear_modes_transitory", - "sim1920_no_shear_modes", - "sim1920_modif_res_no_shear_modes", - "sim3840_modif_res_no_shear_modes", - "sim7680_modif_res_no_shear_modes"] - + directories = [ + "sim960_no_shear_modes", + "sim960_no_shear_modes_transitory", + "sim1920_no_shear_modes", + "sim1920_modif_res_no_shear_modes", + "sim3840_modif_res_no_shear_modes", + "sim7680_modif_res_no_shear_modes", + ] paths_simulations = [] -# directories = [directories[0]] + # directories = [directories[0]] for directory in directories: - paths_simulations += sorted(glob(os.path.join(path_root, directory, "NS2D*"))) + paths_simulations += sorted( + glob(os.path.join(path_root, directory, "NS2D*")) + ) # To compute flow features in all files for path_simulation in paths_simulations: diff --git a/scripts/ns2d.strat/make_table_parameters.py b/scripts/ns2d.strat/make_table_parameters.py index 2c161436..c974c5fd 100644 --- a/scripts/ns2d.strat/make_table_parameters.py +++ b/scripts/ns2d.strat/make_table_parameters.py @@ -30,15 +30,19 @@ path_simulations = sorted(glob(os.path.join(path_root, directory, "NS2D*"))) if MAKE_TABLE: - path_table = ("/home/users/calpelin7m/" + - f"Phd/docs/Manuscript/buoyancy_reynolds_table_n{nx}.tex") - - to_print = ("\\begin{table}[h]\n" - "\\centering \n" - "\\begin{tabular}{cccc} \n" - "\\toprule[1.5pt] \n" + \ - "\\bm{$\\gamma$} & \\bm{$F_h$} & \\bm{$Re_8$} & \\bm{$\\mathcal{R}$} \\\\ \n" - "\\midrule\\ \n") + path_table = ( + "/home/users/calpelin7m/" + + f"Phd/docs/Manuscript/buoyancy_reynolds_table_n{nx}.tex" + ) + + to_print = ( + "\\begin{table}[h]\n" + "\\centering \n" + "\\begin{tabular}{cccc} \n" + "\\toprule[1.5pt] \n" + + "\\bm{$\\gamma$} & \\bm{$F_h$} & \\bm{$Re_8$} & \\bm{$\\mathcal{R}$} \\\\ \n" + "\\midrule\\ \n" + ) for path in path_simulations: @@ -74,31 +78,32 @@ kx = kx[:ikxmax] spectrum1Dkx_EK_ux = spectrum1Dkx_EK_ux[:ikxmax] delta_kx = sim.oper.deltakx - lx = (np.sum(spectrum1Dkx_EK_ux * delta_kx) / - np.sum(kx * spectrum1Dkx_EK_ux * delta_kx)) + lx = np.sum(spectrum1Dkx_EK_ux * delta_kx) / np.sum( + kx * spectrum1Dkx_EK_ux * delta_kx + ) print("lx", lx) # Compute eta_8 - eta_8 = (sim.params.nu_8 ** 3 / epsK_tmean)**(1/22) + eta_8 = (sim.params.nu_8 ** 3 / epsK_tmean) ** (1 / 22) print("eta_8", eta_8) # Compute Re_8 - Re_8 = (lx/eta_8)**(22/3) + Re_8 = (lx / eta_8) ** (22 / 3) print("Re_8", Re_8) # Compute horizontal Froude - F_h = ((epsK_tmean / lx**2)**(1/3)) * (1/sim.params.N) + F_h = ((epsK_tmean / lx ** 2) ** (1 / 3)) * (1 / sim.params.N) print("F_h", F_h) # Reynolds buoyancy 8 - Rb8 = Re_8 * F_h**8 + Rb8 = Re_8 * F_h ** 8 print("Rb8", Rb8) if MAKE_TABLE: - to_print += ("{} & {:.4f} & {:.4e} & {:.4e} \\\\ \n".format( - gamma_table, F_h, Re_8, Rb8)) + to_print += "{} & {:.4f} & {:.4e} & {:.4e} \\\\ \n".format( + gamma_table, F_h, Re_8, Rb8 + ) if MAKE_TABLE: with open(path_table, "w") as f: - to_print += ("\\end{tabular} \n" - "\\end{table}") + to_print += "\\end{tabular} \n" "\\end{table}" f.write(to_print) diff --git a/scripts/ns2d.strat/make_video_spectrum_kykx.py b/scripts/ns2d.strat/make_video_spectrum_kykx.py index 0d7da9d9..f56fea68 100644 --- a/scripts/ns2d.strat/make_video_spectrum_kykx.py +++ b/scripts/ns2d.strat/make_video_spectrum_kykx.py @@ -23,9 +23,10 @@ SAVE = True # Create path with pathlib package. -root_dir = Path( - "/fsnet/project/meige/2015/15DELDUCA/DataSim") / \ - "isotropy_forcing_multidim" +root_dir = ( + Path("/fsnet/project/meige/2015/15DELDUCA/DataSim") + / "isotropy_forcing_multidim" +) list_simulations = [item for item in root_dir.glob("NS2D*")] @@ -51,7 +52,7 @@ skip = 2 tmin = 4 tmax = 400 -scale = "log" # can be "linear" +scale = "log" # can be "linear" # Load object simulation sim = load_state_phys_file(path.as_posix(), merge_missing_params=True) @@ -68,7 +69,6 @@ kx = f["kxE"][...] - ap_fft_spectrum = f["spectrumkykx_ap_fft"] data_ap = ap_fft_spectrum[itmin:itmax:skip, :, :] @@ -77,10 +77,14 @@ # # Create array kz kz = 2 * np.pi * np.fft.fftfreq(poper.ny, poper.Ly / poper.ny) -kz[kz.shape[0]//2] *= -1 +kz[kz.shape[0] // 2] *= -1 kz_modified = np.empty_like(kz) -kz_modified[0:kz_modified.shape[0]//2 - 1] = kz[kz_modified.shape[0]//2 + 1:] -kz_modified[kz_modified.shape[0]//2 - 1:] = kz[0:kz_modified.shape[0]//2 + 1] +kz_modified[0 : kz_modified.shape[0] // 2 - 1] = kz[ + kz_modified.shape[0] // 2 + 1 : +] +kz_modified[kz_modified.shape[0] // 2 - 1 :] = kz[ + 0 : kz_modified.shape[0] // 2 + 1 +] # Create Grid KX, KZ = np.meshgrid(kx, kz_modified) @@ -88,22 +92,26 @@ # Modify data_ap data_ap_plot = data_ap[0, :, :] data_ap_plot_modified = np.empty_like(data_ap_plot) -data_ap_plot_modified[ - 0:kz_modified.shape[0]//2 - 1, :] = data_ap_plot[kz_modified.shape[0]//2 + 1:, :] -data_ap_plot_modified[ - kz_modified.shape[0]//2 - 1:, :] = data_ap_plot[0:kz_modified.shape[0]//2 + 1, :] +data_ap_plot_modified[0 : kz_modified.shape[0] // 2 - 1, :] = data_ap_plot[ + kz_modified.shape[0] // 2 + 1 :, : +] +data_ap_plot_modified[kz_modified.shape[0] // 2 - 1 :, :] = data_ap_plot[ + 0 : kz_modified.shape[0] // 2 + 1, : +] # Modify data_am data_am_plot = data_am[0, :, :] data_am_plot_modified = np.empty_like(data_am_plot) -data_am_plot_modified[ - 0:kz_modified.shape[0]//2 - 1, :] = data_am_plot[kz_modified.shape[0]//2 + 1:, :] -data_am_plot_modified[ - kz_modified.shape[0]//2 - 1:, :] = data_am_plot[0:kz_modified.shape[0]//2 + 1, :] +data_am_plot_modified[0 : kz_modified.shape[0] // 2 - 1, :] = data_am_plot[ + kz_modified.shape[0] // 2 + 1 :, : +] +data_am_plot_modified[kz_modified.shape[0] // 2 - 1 :, :] = data_am_plot[ + 0 : kz_modified.shape[0] // 2 + 1, : +] # Initialize figure fig = plt.figure() -plt.rc('text', usetex=True) +plt.rc("text", usetex=True) ax = fig.add_subplot(221) # ax.set_xlabel(r"$k_x$", fontsize=14) ax.set_ylabel(r"$k_z$", fontsize=14) @@ -129,15 +137,19 @@ forcing_rate = pforcing.forcing_rate l_f = 2 * np.pi / (pforcing.nkmax_forcing * sim.oper.deltaky) # energy_f = 100 * (forcing_rate ** 6 * l_f**2)**(1/7) -energy_f = ((forcing_rate ** 2) * (l_f**10))**(1/7) +energy_f = ((forcing_rate ** 2) * (l_f ** 10)) ** (1 / 7) # Compute energy energies_ap = np.empty_like(times) energies_am = np.empty_like(times) for it, time in enumerate(times): - energies_ap[it] = (sim.oper.deltakx * sim.oper.deltaky) * data_ap[it].sum() / energy_f - energies_am[it] = (sim.oper.deltakx * sim.oper.deltaky) * data_am[it].sum() / energy_f + energies_ap[it] = ( + (sim.oper.deltakx * sim.oper.deltaky) * data_ap[it].sum() / energy_f + ) + energies_am[it] = ( + (sim.oper.deltakx * sim.oper.deltaky) * data_am[it].sum() / energy_f + ) # Set axes limit xlim = 200 @@ -166,26 +178,32 @@ # Plot first figure -ax.plot(kx, 1e-0 * sim.params.N * (kx / eps)**(1/3), color="white") -ax.plot(kx, -1e-0 * sim.params.N * (kx / eps)**(1/3), color="white") +ax.plot(kx, 1e-0 * sim.params.N * (kx / eps) ** (1 / 3), color="white") +ax.plot(kx, -1e-0 * sim.params.N * (kx / eps) ** (1 / 3), color="white") ax.text(kx[ikx_text], kz[ikz_text], r"\hat{a}_+", color="white", fontsize=15) -ax2.plot(kx, 1e-0 * sim.params.N * (kx / eps)**(1/3), color="white") -ax2.plot(kx, -1e-0 * sim.params.N * (kx / eps)**(1/3), color="white") +ax2.plot(kx, 1e-0 * sim.params.N * (kx / eps) ** (1 / 3), color="white") +ax2.plot(kx, -1e-0 * sim.params.N * (kx / eps) ** (1 / 3), color="white") ax2.text(kx[ikx_text], kz[ikz_text], r"\hat{a}_-", color="white", fontsize=15) if scale == "linear": - _im = ax.pcolormesh(KX, KZ, data_ap_plot_modified, vmin=vmin, vmax=vmax) - _im2 = ax2.pcolormesh(KX, KZ, data_am_plot_modified, vmin=vmin, vmax=vmax) + _im = ax.pcolormesh(KX, KZ, data_ap_plot_modified, vmin=vmin, vmax=vmax) + _im2 = ax2.pcolormesh(KX, KZ, data_am_plot_modified, vmin=vmin, vmax=vmax) elif scale == "log": - _im = ax.pcolormesh(KX, KZ, np.log10(data_ap_plot_modified), vmin=vmin, vmax=vmax) - _im2 = ax2.pcolormesh(KX, KZ, np.log10(data_am_plot_modified), vmin=vmin, vmax=vmax) + _im = ax.pcolormesh( + KX, KZ, np.log10(data_ap_plot_modified), vmin=vmin, vmax=vmax + ) + _im2 = ax2.pcolormesh( + KX, KZ, np.log10(data_am_plot_modified), vmin=vmin, vmax=vmax + ) ax1.plot(times, energies_ap, color="grey", alpha=0.4) ax1.plot(times, energies_am, color="grey", alpha=0.4) _im_inset = ax1.plot(times[0], energies_ap[0], color="red", label=r"$\hat{a}_+$") -_im_inset3 = ax1.plot(times[0], energies_am[0], color="blue", label=r"$\hat{a}_-$") +_im_inset3 = ax1.plot( + times[0], energies_am[0], color="blue", label=r"$\hat{a}_-$" +) ax1.legend(fontsize=14) cbar_ax = fig.add_axes([0.38, 0.15, 0.01, 0.7]) colorbar = fig.colorbar(_im, cax=cbar_ax, format="%.1f") @@ -196,33 +214,35 @@ def _update(frame): data_ap_plot = data_ap[frame, :, :] data_ap_plot_modified = np.empty_like(data_ap_plot) - data_ap_plot_modified[0:kz_modified.shape[0]//2 - 1, :] = data_ap_plot[ - kz_modified.shape[0]//2 + 1:, :] - data_ap_plot_modified[kz_modified.shape[0]//2 - 1:, :] = data_ap_plot[ - 0:kz_modified.shape[0]//2 + 1, :] - + data_ap_plot_modified[0 : kz_modified.shape[0] // 2 - 1, :] = data_ap_plot[ + kz_modified.shape[0] // 2 + 1 :, : + ] + data_ap_plot_modified[kz_modified.shape[0] // 2 - 1 :, :] = data_ap_plot[ + 0 : kz_modified.shape[0] // 2 + 1, : + ] data_am_plot = data_am[frame, :, :] data_am_plot_modified = np.empty_like(data_am_plot) - data_am_plot_modified[0:kz_modified.shape[0]//2 - 1, :] = data_am_plot[ - kz_modified.shape[0]//2 + 1:, :] - data_am_plot_modified[kz_modified.shape[0]//2 - 1:, :] = data_am_plot[ - 0:kz_modified.shape[0]//2 + 1, :] + data_am_plot_modified[0 : kz_modified.shape[0] // 2 - 1, :] = data_am_plot[ + kz_modified.shape[0] // 2 + 1 :, : + ] + data_am_plot_modified[kz_modified.shape[0] // 2 - 1 :, :] = data_am_plot[ + 0 : kz_modified.shape[0] // 2 + 1, : + ] # Trick: https://stackoverflow.com/questions/29009743/using-set-array-with-pyplot-pcolormesh-ruins-figure if scale == "linear": - _im.set_clim(vmin=0, vmax= 0.001 * data_ap_plot_modified.max()) + _im.set_clim(vmin=0, vmax=0.001 * data_ap_plot_modified.max()) _im.set_array(data_ap_plot_modified[:-1, :-1].flatten()) - _im2.set_clim(vmin=0, vmax= 0.001 * data_ap_plot_modified.max()) + _im2.set_clim(vmin=0, vmax=0.001 * data_ap_plot_modified.max()) _im2.set_array(data_am_plot_modified[:-1, :-1].flatten()) elif scale == "log": - _im.set_clim(vmin=-5, vmax= np.log10(data_ap_plot_modified.max()) - 3) + _im.set_clim(vmin=-5, vmax=np.log10(data_ap_plot_modified.max()) - 3) _im.set_array(np.log10(data_ap_plot_modified[:-1, :-1].flatten())) - - _im2.set_clim(vmin=-5, vmax= np.log10(data_ap_plot_modified.max()) - 3) + _im2.set_clim(vmin=-5, vmax=np.log10(data_ap_plot_modified.max()) - 3) _im2.set_array(np.log10(data_am_plot_modified[:-1, :-1].flatten())) _im_inset[0].set_data(times[0:frame], energies_ap[0:frame]) @@ -230,12 +250,16 @@ def _update(frame): ax.set_title(r"$t = {:.0f} \tau_{{af}}$".format(times[frame]), fontsize=14) + ani = animation.FuncAnimation( - fig, _update, len(times), interval=1000, repeat=True) + fig, _update, len(times), interval=1000, repeat=True +) if SAVE: - ani.save("/home/users/calpelin7m/Phd/Movies_spectrakzkx/spectrumkykx_{}_kznegative_{}.mp4".format( - sim.params.forcing.key_forced, - bool(sim.params.forcing.tcrandom_anisotropic.kz_negative_enable)), - writer="ffmpeg" + ani.save( + "/home/users/calpelin7m/Phd/Movies_spectrakzkx/spectrumkykx_{}_kznegative_{}.mp4".format( + sim.params.forcing.key_forced, + bool(sim.params.forcing.tcrandom_anisotropic.kz_negative_enable), + ), + writer="ffmpeg", ) diff --git a/scripts/ns2d.strat/ns2dstrat_lmode.py b/scripts/ns2d.strat/ns2dstrat_lmode.py index d90bd767..90425f8b 100644 --- a/scripts/ns2d.strat/ns2dstrat_lmode.py +++ b/scripts/ns2d.strat/ns2dstrat_lmode.py @@ -9,14 +9,17 @@ from math import pi from fluidsim.solvers.ns2d.strat.solver import Simul -def make_parameters_simulation(gamma, F, sigma, nu_8, t_end=10, NO_SHEAR_MODES=False): + +def make_parameters_simulation( + gamma, F, sigma, nu_8, t_end=10, NO_SHEAR_MODES=False +): ## Operator parameters - anisotropy_domain = 4 # anisotropy_domain = nx / nz + anisotropy_domain = 4 # anisotropy_domain = nx / nz nx = 240 nz = nx // anisotropy_domain Lx = 2 * pi - Lz = Lx * (nz / nx) # deltax = deltay + Lz = Lx * (nz / nx) # deltax = deltay coef_dealiasing = 0.6666 @@ -29,7 +32,7 @@ def make_parameters_simulation(gamma, F, sigma, nu_8, t_end=10, NO_SHEAR_MODES=F forcing_enable = True nkmax_forcing = 8 nkmin_forcing = 4 - tau_af = 1 # Forcing time equal to 1 + tau_af = 1 # Forcing time equal to 1 ###################### ####################### @@ -47,7 +50,7 @@ def make_parameters_simulation(gamma, F, sigma, nu_8, t_end=10, NO_SHEAR_MODES=F # Forcing parameters params.forcing.enable = forcing_enable - params.forcing.type = 'tcrandom_anisotropic' + params.forcing.type = "tcrandom_anisotropic" params.forcing.key_forced = "ap_fft" params.forcing.nkmax_forcing = nkmax_forcing params.forcing.nkmin_forcing = nkmin_forcing @@ -56,14 +59,16 @@ def make_parameters_simulation(gamma, F, sigma, nu_8, t_end=10, NO_SHEAR_MODES=F # Compute other parameters k_f = ((nkmax_forcing + nkmin_forcing) / 2) * max(2 * pi / Lx, 2 * pi / Lz) - forcing_rate = (1 / tau_af**7) * ((2 * pi) / k_f)**2 + forcing_rate = (1 / tau_af ** 7) * ((2 * pi) / k_f) ** 2 omega_af = 2 * pi / tau_af params.N = (gamma / F) * omega_af params.nu_8 = nu_8 # Continuation on forcing... params.forcing.forcing_rate = forcing_rate - params.forcing.tcrandom.time_correlation = sigma * (pi / (params.N * F)) # time_correlation = wave period + params.forcing.tcrandom.time_correlation = sigma * ( + pi / (params.N * F) + ) # time_correlation = wave period # Time stepping parameters params.time_stepping.USE_CFL = USE_CFL @@ -77,6 +82,7 @@ def make_parameters_simulation(gamma, F, sigma, nu_8, t_end=10, NO_SHEAR_MODES=F return params + def modify_parameters(params): # Output parameters @@ -86,16 +92,17 @@ def modify_parameters(params): params.output.periods_save.spect_energy_budg = 1e-1 params.output.periods_save.spectra = 1e-1 + if __name__ == "__main__": ##### PARAMETERS ##### ###################### - gamma = 0.2 # gamma = omega_l / omega_af - F = np.sin(pi / 4) # F = omega_l / N - sigma = 1 # sigma = omega_l / (pi * f_cf); f_cf freq time correlation forcing in s-1 + gamma = 0.2 # gamma = omega_l / omega_af + F = np.sin(pi / 4) # F = omega_l / N + sigma = 1 # sigma = omega_l / (pi * f_cf); f_cf freq time correlation forcing in s-1 nu_8 = 1e-15 - params = make_parameters_simulation(gamma, F, sigma, nu_8, t_end=50.) + params = make_parameters_simulation(gamma, F, sigma, nu_8, t_end=50.0) # Start simulation sim = Simul(params) @@ -105,19 +112,19 @@ def modify_parameters(params): # Normalize initialization if sim.params.oper.nx != sim.params.oper.ny: KX = sim.oper.KX - cond = KX == 0. + cond = KX == 0.0 - ux_fft = sim.state.get_var('ux_fft') - uy_fft = sim.state.get_var('uy_fft') - b_fft = sim.state.get_var('b_fft') + ux_fft = sim.state.get_var("ux_fft") + uy_fft = sim.state.get_var("uy_fft") + b_fft = sim.state.get_var("b_fft") - ux_fft[cond] = 0. - uy_fft[cond] = 0. - b_fft[cond] = 0. + ux_fft[cond] = 0.0 + uy_fft[cond] = 0.0 + b_fft[cond] = 0.0 # Compute energy after ux_fft[kx=0] uy_fft[kx=0] b_fft[kx=0] - ek_fft = (np.abs(ux_fft)**2 + np.abs(uy_fft)**2)/2 - ea_fft = ((np.abs(b_fft)/params.N)**2)/2 + ek_fft = (np.abs(ux_fft) ** 2 + np.abs(uy_fft) ** 2) / 2 + ea_fft = ((np.abs(b_fft) / params.N) ** 2) / 2 e_fft = ek_fft + ea_fft energy_before_norm = sim.output.sum_wavenumbers(e_fft) @@ -128,8 +135,10 @@ def modify_parameters(params): nkmax_forcing = params.forcing.nkmax_forcing nkmin_forcing = params.forcing.nkmin_forcing - k_f = ((nkmax_forcing + nkmin_forcing) / 2) * max(2 * pi / Lx, 2 * pi / Lz) - energy_f = params.forcing.forcing_rate**(2/7) * (2 * pi / k_f)**7 + k_f = ((nkmax_forcing + nkmin_forcing) / 2) * max( + 2 * pi / Lx, 2 * pi / Lz + ) + energy_f = params.forcing.forcing_rate ** (2 / 7) * (2 * pi / k_f) ** 7 coef = np.sqrt(1e-4 * energy_f / energy_before_norm) diff --git a/scripts/ns2d.strat/plot_buoyancy_reynolds_anisotropy.py b/scripts/ns2d.strat/plot_buoyancy_reynolds_anisotropy.py index 161ffbe6..6383f184 100644 --- a/scripts/ns2d.strat/plot_buoyancy_reynolds_anisotropy.py +++ b/scripts/ns2d.strat/plot_buoyancy_reynolds_anisotropy.py @@ -17,12 +17,14 @@ from fluiddyn.output.rcparams import set_rcparams path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" -directories = ["sim960_no_shear_modes", - "sim960_no_shear_modes_transitory", - "sim1920_no_shear_modes", - "sim1920_modif_res_no_shear_modes", - "sim3840_modif_res_no_shear_modes", - "sim7680_modif_res_no_shear_modes"] +directories = [ + "sim960_no_shear_modes", + "sim960_no_shear_modes_transitory", + "sim1920_no_shear_modes", + "sim1920_modif_res_no_shear_modes", + "sim3840_modif_res_no_shear_modes", + "sim7680_modif_res_no_shear_modes", +] paths_simulations = [] @@ -77,35 +79,75 @@ cax = divider.append_axes("right", size="5%", pad=0.05) cax.tick_params(labelsize=14) fig.colorbar(scatter, cax=cax) -ax.text(1e7, 0.013, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=12) +ax.text( + 1e7, + 0.013, + r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", + fontsize=12, +) # Legend... -blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None', - markersize=8, label=r'$n_x = 960$') -red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None', - markersize=8, label=r'$n_x = 1920$') -purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 3840$') -diamond = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 7680$') - -ax.legend(handles=[blue_star, red_square, purple_triangle, diamond], - loc="upper center", - bbox_to_anchor=(0.5,1.1), - borderaxespad=0., - ncol=len(markers), - handletextpad=0.1, - fontsize=14) +blue_star = mlines.Line2D( + [], + [], + color="red", + marker="o", + linestyle="None", + markersize=8, + label=r"$n_x = 960$", +) +red_square = mlines.Line2D( + [], + [], + color="red", + marker="s", + linestyle="None", + markersize=8, + label=r"$n_x = 1920$", +) +purple_triangle = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 3840$", +) +diamond = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 7680$", +) + +ax.legend( + handles=[blue_star, red_square, purple_triangle, diamond], + loc="upper center", + bbox_to_anchor=(0.5, 1.1), + borderaxespad=0.0, + ncol=len(markers), + handletextpad=0.1, + fontsize=14, +) # Add text... -ax.text(1e-6, 1e-1, - r"$A = \frac{\log_{10} \left( 2E_{k,x} / E_k \right)}{\log_{10}(2)}$", - fontsize=20) +ax.text( + 1e-6, + 1e-1, + r"$A = \frac{\log_{10} \left( 2E_{k,x} / E_k \right)}{\log_{10}(2)}$", + fontsize=20, +) SAVE = True if SAVE: - path_save = "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + \ - "buoyancy_reynolds_anisotropy.png" + path_save = ( + "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + + "buoyancy_reynolds_anisotropy.png" + ) fig.savefig(path_save, format="png", bbox_inches="tight") diff --git a/scripts/ns2d.strat/plot_buoyancy_reynolds_dissipation.py b/scripts/ns2d.strat/plot_buoyancy_reynolds_dissipation.py index 2930b5b2..fa41adc2 100644 --- a/scripts/ns2d.strat/plot_buoyancy_reynolds_dissipation.py +++ b/scripts/ns2d.strat/plot_buoyancy_reynolds_dissipation.py @@ -17,12 +17,14 @@ from fluiddyn.output.rcparams import set_rcparams path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" -directories = ["sim960_no_shear_modes", - "sim960_no_shear_modes_transitory", - "sim1920_no_shear_modes", - "sim1920_modif_res_no_shear_modes", - "sim3840_modif_res_no_shear_modes", - "sim7680_modif_res_no_shear_modes"] +directories = [ + "sim960_no_shear_modes", + "sim960_no_shear_modes_transitory", + "sim1920_no_shear_modes", + "sim1920_modif_res_no_shear_modes", + "sim3840_modif_res_no_shear_modes", + "sim7680_modif_res_no_shear_modes", +] paths_simulations = [] @@ -80,22 +82,52 @@ # ax.text(1e7, 0.013, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=12) # Legend... -blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None', - markersize=8, label=r'$n_x = 960$') -red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None', - markersize=8, label=r'$n_x = 1920$') -purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 3840$') -diamond = mlines.Line2D([], [], color='red', marker='d', linestyle='None', - markersize=8, label=r'$n_x = 7680$') - -ax.legend(handles=[blue_star, red_square, purple_triangle, diamond], - loc="upper center", - bbox_to_anchor=(0.5,1.1), - borderaxespad=0., - ncol=len(markers), - handletextpad=0.1, - fontsize=14) +blue_star = mlines.Line2D( + [], + [], + color="red", + marker="o", + linestyle="None", + markersize=8, + label=r"$n_x = 960$", +) +red_square = mlines.Line2D( + [], + [], + color="red", + marker="s", + linestyle="None", + markersize=8, + label=r"$n_x = 1920$", +) +purple_triangle = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 3840$", +) +diamond = mlines.Line2D( + [], + [], + color="red", + marker="d", + linestyle="None", + markersize=8, + label=r"$n_x = 7680$", +) + +ax.legend( + handles=[blue_star, red_square, purple_triangle, diamond], + loc="upper center", + bbox_to_anchor=(0.5, 1.1), + borderaxespad=0.0, + ncol=len(markers), + handletextpad=0.1, + fontsize=14, +) # Add text... # ax.text(1e-6, 1e-1, @@ -104,8 +136,10 @@ SAVE = True if SAVE: - path_save = "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + \ - "buoyancy_reynolds_dissipation.png" + path_save = ( + "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + + "buoyancy_reynolds_dissipation.png" + ) fig.savefig(path_save, format="png", bbox_inches="tight") diff --git a/scripts/ns2d.strat/plot_buoyancy_reynolds_froude.py b/scripts/ns2d.strat/plot_buoyancy_reynolds_froude.py index ef65d1b4..35248dd6 100644 --- a/scripts/ns2d.strat/plot_buoyancy_reynolds_froude.py +++ b/scripts/ns2d.strat/plot_buoyancy_reynolds_froude.py @@ -14,11 +14,13 @@ from fluiddyn.output.rcparams import set_rcparams path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" -directories = ["sim960_no_shear_modes", - "sim960_no_shear_modes_transitory", - "sim1920_no_shear_modes", - "sim1920_modif_res_no_shear_modes", - "sim7680_modif_res_no_shear_modes"] +directories = [ + "sim960_no_shear_modes", + "sim960_no_shear_modes_transitory", + "sim1920_no_shear_modes", + "sim1920_modif_res_no_shear_modes", + "sim7680_modif_res_no_shear_modes", +] paths_simulations = [] @@ -76,27 +78,59 @@ # ax.text(1e7, 0.013, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=12) # Legend... -blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None', - markersize=8, label=r'$n_x = 960$') -red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None', - markersize=8, label=r'$n_x = 1920$') -purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 3840$') -diamond = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 7680$') - -ax.legend(handles=[blue_star, red_square, purple_triangle, diamond], - loc="upper center", - bbox_to_anchor=(0.5,1.1), - borderaxespad=0., - ncol=len(markers), - handletextpad=0.1, - fontsize=14) +blue_star = mlines.Line2D( + [], + [], + color="red", + marker="o", + linestyle="None", + markersize=8, + label=r"$n_x = 960$", +) +red_square = mlines.Line2D( + [], + [], + color="red", + marker="s", + linestyle="None", + markersize=8, + label=r"$n_x = 1920$", +) +purple_triangle = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 3840$", +) +diamond = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 7680$", +) + +ax.legend( + handles=[blue_star, red_square, purple_triangle, diamond], + loc="upper center", + bbox_to_anchor=(0.5, 1.1), + borderaxespad=0.0, + ncol=len(markers), + handletextpad=0.1, + fontsize=14, +) SAVE = True if SAVE: - path_save = "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + \ - "buoyancy_reynolds_froude.png" + path_save = ( + "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + + "buoyancy_reynolds_froude.png" + ) fig.savefig(path_save, format="png", bbox_inches="tight") diff --git a/scripts/ns2d.strat/plot_froude_anisotropy.py b/scripts/ns2d.strat/plot_froude_anisotropy.py index 04e60cb5..b3ac106e 100644 --- a/scripts/ns2d.strat/plot_froude_anisotropy.py +++ b/scripts/ns2d.strat/plot_froude_anisotropy.py @@ -17,12 +17,14 @@ from fluiddyn.output.rcparams import set_rcparams path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" -directories = ["sim960_no_shear_modes", - "sim960_no_shear_modes_transitory", - "sim1920_no_shear_modes", - "sim1920_modif_res_no_shear_modes", - "sim3840_modif_res_no_shear_modes", - "sim7680_modif_res_no_shear_modes"] +directories = [ + "sim960_no_shear_modes", + "sim960_no_shear_modes_transitory", + "sim1920_no_shear_modes", + "sim1920_modif_res_no_shear_modes", + "sim3840_modif_res_no_shear_modes", + "sim7680_modif_res_no_shear_modes", +] paths_simulations = [] @@ -76,34 +78,74 @@ cax = divider.append_axes("right", size="5%", pad=0.05) cax.tick_params(labelsize=14) fig.colorbar(scatter, cax=cax) -ax.text(1.1, 0.013, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=12) +ax.text( + 1.1, + 0.013, + r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", + fontsize=12, +) # Legend... -blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None', - markersize=8, label=r'$n_x = 960$') -red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None', - markersize=8, label=r'$n_x = 1920$') -purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 3840$') -diamond = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 7680$') -ax.legend(handles=[blue_star, red_square, purple_triangle, diamond], - loc="upper center", - bbox_to_anchor=(0.5,1.1), - borderaxespad=0., - ncol=len(markers), - handletextpad=0.1, - fontsize=14) +blue_star = mlines.Line2D( + [], + [], + color="red", + marker="o", + linestyle="None", + markersize=8, + label=r"$n_x = 960$", +) +red_square = mlines.Line2D( + [], + [], + color="red", + marker="s", + linestyle="None", + markersize=8, + label=r"$n_x = 1920$", +) +purple_triangle = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 3840$", +) +diamond = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 7680$", +) +ax.legend( + handles=[blue_star, red_square, purple_triangle, diamond], + loc="upper center", + bbox_to_anchor=(0.5, 1.1), + borderaxespad=0.0, + ncol=len(markers), + handletextpad=0.1, + fontsize=14, +) # Add text... -ax.text(1e-2, 1e-1, - r"$A = \frac{\log_{10} \left( 2E_{k,x} / E_k \right)}{\log_{10}(2)}$", - fontsize=20) +ax.text( + 1e-2, + 1e-1, + r"$A = \frac{\log_{10} \left( 2E_{k,x} / E_k \right)}{\log_{10}(2)}$", + fontsize=20, +) SAVE = True if SAVE: - path_save = "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + \ - "froude_anisotropy.png" + path_save = ( + "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + + "froude_anisotropy.png" + ) fig.savefig(path_save, format="png", bbox_inches="tight") diff --git a/scripts/ns2d.strat/plot_froude_dissipation.py b/scripts/ns2d.strat/plot_froude_dissipation.py index 6845b9a8..5538dc4b 100644 --- a/scripts/ns2d.strat/plot_froude_dissipation.py +++ b/scripts/ns2d.strat/plot_froude_dissipation.py @@ -17,12 +17,14 @@ from fluiddyn.output.rcparams import set_rcparams path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" -directories = ["sim960_no_shear_modes", - "sim960_no_shear_modes_transitory", - "sim1920_no_shear_modes", - "sim1920_modif_res_no_shear_modes", - "sim3840_modif_res_no_shear_modes", - "sim7680_modif_res_no_shear_modes"] +directories = [ + "sim960_no_shear_modes", + "sim960_no_shear_modes_transitory", + "sim1920_no_shear_modes", + "sim1920_modif_res_no_shear_modes", + "sim3840_modif_res_no_shear_modes", + "sim7680_modif_res_no_shear_modes", +] paths_simulations = [] @@ -80,22 +82,52 @@ # ax.text(1e7, 0.013, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=12) # Legend... -blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None', - markersize=8, label=r'$n_x = 960$') -red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None', - markersize=8, label=r'$n_x = 1920$') -purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 3840$') -diamond = mlines.Line2D([], [], color='red', marker='d', linestyle='None', - markersize=8, label=r'$n_x = 7680$') - -ax.legend(handles=[blue_star, red_square, purple_triangle, diamond], - loc="upper center", - bbox_to_anchor=(0.5,1.1), - borderaxespad=0., - ncol=len(markers), - handletextpad=0.1, - fontsize=14) +blue_star = mlines.Line2D( + [], + [], + color="red", + marker="o", + linestyle="None", + markersize=8, + label=r"$n_x = 960$", +) +red_square = mlines.Line2D( + [], + [], + color="red", + marker="s", + linestyle="None", + markersize=8, + label=r"$n_x = 1920$", +) +purple_triangle = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 3840$", +) +diamond = mlines.Line2D( + [], + [], + color="red", + marker="d", + linestyle="None", + markersize=8, + label=r"$n_x = 7680$", +) + +ax.legend( + handles=[blue_star, red_square, purple_triangle, diamond], + loc="upper center", + bbox_to_anchor=(0.5, 1.1), + borderaxespad=0.0, + ncol=len(markers), + handletextpad=0.1, + fontsize=14, +) # Add text... # ax.text(1e-6, 1e-1, @@ -104,8 +136,10 @@ SAVE = True if SAVE: - path_save = "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + \ - "froude_dissipation.png" + path_save = ( + "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + + "froude_dissipation.png" + ) fig.savefig(path_save, format="png", bbox_inches="tight") diff --git a/scripts/ns2d.strat/plot_length_scales.py b/scripts/ns2d.strat/plot_length_scales.py index 14c3dea0..f7c1d5aa 100644 --- a/scripts/ns2d.strat/plot_length_scales.py +++ b/scripts/ns2d.strat/plot_length_scales.py @@ -10,6 +10,7 @@ from glob import glob from mpl_toolkits.axes_grid1 import make_axes_locatable + # from compute_anisotropy import compute_anisotropy # from compute_ratio_dissipation import compute_ratio_dissipation # from compute_reynolds_froude import compute_buoyancy_reynolds @@ -19,20 +20,26 @@ from fluidsim import load_params_simul from flow_features import get_features_from_sim, _get_resolution_from_dir + def _get_resolution_from_dir(path_simulation): return path_simulation.split("NS2D.strat_")[1].split("x")[0] + + def _get_gamma_str_from_path(path_simulation): return path_simulation.split("_gamma")[1].split("_")[0] + SAVE = False # Create path simulations path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" -directories = ["sim960_no_shear_modes", - "sim960_no_shear_modes_transitory", - "sim1920_no_shear_modes", - "sim1920_modif_res_no_shear_modes", - "sim3840_modif_res_no_shear_modes"] +directories = [ + "sim960_no_shear_modes", + "sim960_no_shear_modes_transitory", + "sim1920_no_shear_modes", + "sim1920_modif_res_no_shear_modes", + "sim3840_modif_res_no_shear_modes", +] # directories = ["sim960_no_shear_modes", # "sim1920_no_shear_modes", @@ -84,7 +91,7 @@ def _get_gamma_str_from_path(path_simulation): for path_file in path_phys_files[-10:]: with h5py.File(path_file, "r") as f: ux = f["state_phys"]["ux"][...] - ux_rms.append(np.sqrt(np.mean(ux**2))) + ux_rms.append(np.sqrt(np.mean(ux ** 2))) # Load parameters params = load_params_simul(path) @@ -108,8 +115,19 @@ def _get_gamma_str_from_path(path_simulation): print("Re_8", Re_8) print("R_b", R_b) - for _f, _r, _a, _d, _m, _l in zip(froudes, reynoldsb, anisotropies, dissipations, markers, lzs_billant): - scatter = ax.scatter(_r, _l, s=100 * (1**1), c="k", vmin=0, vmax=1, marker=_m, alpha=0.7) + for _f, _r, _a, _d, _m, _l in zip( + froudes, reynoldsb, anisotropies, dissipations, markers, lzs_billant + ): + scatter = ax.scatter( + _r, + _l, + s=100 * (1 ** 1), + c="k", + vmin=0, + vmax=1, + marker=_m, + alpha=0.7, + ) # for _f, _r, _a, _d, _m, _l in zip(froudes, reynoldsb, anisotropies, dissipations, markers, lzs_billant): # scatter = ax.scatter(_r, _l, s=2000 * (_a**1) + 20, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7) @@ -132,25 +150,51 @@ def _get_gamma_str_from_path(path_simulation): import matplotlib.lines as mlines import matplotlib.pyplot as plt -blue_star = mlines.Line2D([], [], color='k', marker='o', linestyle='None', - markersize=8, label=r'$n_x = 960$') -red_square = mlines.Line2D([], [], color='k', marker='s', linestyle='None', - markersize=8, label=r'$n_x = 1920$') -purple_triangle = mlines.Line2D([], [], color='k', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 3840$') - -ax.legend(handles=[blue_star, red_square, purple_triangle], - loc="upper center", - bbox_to_anchor=(0.5,1.1), - borderaxespad=0., - ncol=len(markers), - handletextpad=0.1, - fontsize=14) +blue_star = mlines.Line2D( + [], + [], + color="k", + marker="o", + linestyle="None", + markersize=8, + label=r"$n_x = 960$", +) +red_square = mlines.Line2D( + [], + [], + color="k", + marker="s", + linestyle="None", + markersize=8, + label=r"$n_x = 1920$", +) +purple_triangle = mlines.Line2D( + [], + [], + color="k", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 3840$", +) + +ax.legend( + handles=[blue_star, red_square, purple_triangle], + loc="upper center", + bbox_to_anchor=(0.5, 1.1), + borderaxespad=0.0, + ncol=len(markers), + handletextpad=0.1, + fontsize=14, +) fig.tight_layout(pad=0.4) if SAVE: path_save = "/home/users/calpelin7m/Phd/docs/EFMC18/figures" - fig.savefig(path_save + "/vertical_length_billant.eps", format="eps", - bbox_inches="tight") + fig.savefig( + path_save + "/vertical_length_billant.eps", + format="eps", + bbox_inches="tight", + ) plt.show() diff --git a/scripts/ns2d.strat/plot_reynolds_froude.py b/scripts/ns2d.strat/plot_reynolds_froude.py index c1522dfd..1ae8b3e3 100644 --- a/scripts/ns2d.strat/plot_reynolds_froude.py +++ b/scripts/ns2d.strat/plot_reynolds_froude.py @@ -17,12 +17,15 @@ from compute_reynolds_froude import compute_buoyancy_reynolds from fluiddyn.output.rcparams import set_rcparams + def _get_resolution_from_dir(path_simulation): return path_simulation.split("NS2D.strat_")[1].split("x")[0] + def _get_gamma_str_from_path(path_simulation): return path_simulation.split("_gamma")[1].split("_")[0] + SAVE = False # Create path simulations @@ -54,7 +57,12 @@ def _get_gamma_str_from_path(path_simulation): ax.set_ylabel(r"$\mathcal{R}_8$", fontsize=18) ax.set_xscale("log") ax.set_yscale("log") -ax.text(0.6, 2e-12, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=16) +ax.text( + 0.6, + 2e-12, + r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", + fontsize=16, +) ax.tick_params(axis="x", labelsize=18) ax.tick_params(axis="y", labelsize=18) ax.set_xlim([0.003, 1]) @@ -91,12 +99,27 @@ def _get_gamma_str_from_path(path_simulation): # for _f, _r, _a, _d, _m in zip(froudes, reynoldsb, anisotropies, dissipations, markers): # scatter = ax.scatter(_f, _r, s=2000 * (_a**1) + 20, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7) -for _f, _r, _a, _d, _m in zip(froudes, reynoldsb, anisotropies, dissipations, markers): - scatter = ax.scatter(_f, _r, s=2000 * (_a**1) + 20, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7) - -ax.scatter(0.1, 1e-5, marker="o", s=2000 * (np.log10(2*0.5)**1) + 20, color="red") -ax.scatter(0.1, 1e-6, marker="o", s=2000 * (np.log10(2*0.75)**1) + 20, color="red") -ax.scatter(0.1, 2e-8, marker="o", s=2000 * (np.log10(2)**1) + 20, color="red") +for _f, _r, _a, _d, _m in zip( + froudes, reynoldsb, anisotropies, dissipations, markers +): + scatter = ax.scatter( + _f, + _r, + s=2000 * (_a ** 1) + 20, + c=_d, + vmin=0, + vmax=1, + marker=_m, + alpha=0.7, + ) + +ax.scatter( + 0.1, 1e-5, marker="o", s=2000 * (np.log10(2 * 0.5) ** 1) + 20, color="red" +) +ax.scatter( + 0.1, 1e-6, marker="o", s=2000 * (np.log10(2 * 0.75) ** 1) + 20, color="red" +) +ax.scatter(0.1, 2e-8, marker="o", s=2000 * (np.log10(2) ** 1) + 20, color="red") ax.text(0.16, 8e-6, r"$anisotropy = 0$", color="red", fontsize=14) ax.text(0.16, 5e-7, r"$anisotropy = 0.5$", color="red", fontsize=14) ax.text(0.16, 2e-8, r"$anisotropy = 1$", color="red", fontsize=14) @@ -109,25 +132,49 @@ def _get_gamma_str_from_path(path_simulation): import matplotlib.lines as mlines import matplotlib.pyplot as plt -blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None', - markersize=8, label=r'$n_x = 960$') -red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None', - markersize=8, label=r'$n_x = 1920$') -purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 3840$') - -ax.legend(handles=[blue_star, red_square, purple_triangle], - loc="upper center", - bbox_to_anchor=(0.5,1.1), - borderaxespad=0., - ncol=len(markers), - handletextpad=0.1, - fontsize=14) +blue_star = mlines.Line2D( + [], + [], + color="red", + marker="o", + linestyle="None", + markersize=8, + label=r"$n_x = 960$", +) +red_square = mlines.Line2D( + [], + [], + color="red", + marker="s", + linestyle="None", + markersize=8, + label=r"$n_x = 1920$", +) +purple_triangle = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 3840$", +) + +ax.legend( + handles=[blue_star, red_square, purple_triangle], + loc="upper center", + bbox_to_anchor=(0.5, 1.1), + borderaxespad=0.0, + ncol=len(markers), + handletextpad=0.1, + fontsize=14, +) fig.tight_layout(pad=0.4) if SAVE: path_save = "/home/users/calpelin7m/Phd/docs/EFMC18/figures" - fig.savefig(path_save + "/reynoldsb_froude.eps", format="eps", - bbox_inches="tight") + fig.savefig( + path_save + "/reynoldsb_froude.eps", format="eps", bbox_inches="tight" + ) plt.show() diff --git a/scripts/ns2d.strat/plot_vertical_scale.py b/scripts/ns2d.strat/plot_vertical_scale.py index f9439edb..6831336d 100644 --- a/scripts/ns2d.strat/plot_vertical_scale.py +++ b/scripts/ns2d.strat/plot_vertical_scale.py @@ -20,12 +20,14 @@ from fluiddyn.output.rcparams import set_rcparams path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" -directories = ["sim960_no_shear_modes", - "sim960_no_shear_modes_transitory", - "sim1920_no_shear_modes", - "sim1920_modif_res_no_shear_modes", - "sim3840_modif_res_no_shear_modes", - "sim7680_modif_res_no_shear_modes"] +directories = [ + "sim960_no_shear_modes", + "sim960_no_shear_modes_transitory", + "sim1920_no_shear_modes", + "sim1920_modif_res_no_shear_modes", + "sim3840_modif_res_no_shear_modes", + "sim7680_modif_res_no_shear_modes", +] paths_simulations = [] @@ -68,7 +70,7 @@ for path_file in path_phys_files[-10:-1]: with h5py.File(path_file, "r") as f: ux = f["state_phys"]["ux"][...] - ux_rms.append(np.sqrt(np.mean(ux**2))) + ux_rms.append(np.sqrt(np.mean(ux ** 2))) # Load parameters params = load_params_simul(path) @@ -86,8 +88,12 @@ ax.set_ylim([0, 1.2e1]) # Plot... -for _f, _r, _a, _d, _m, _l in zip(froudes, reynoldsb, anisotropies, dissipations, markers, lzs_billant): - scatter = ax.scatter(_r, _l, s=100, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7) +for _f, _r, _a, _d, _m, _l in zip( + froudes, reynoldsb, anisotropies, dissipations, markers, lzs_billant +): + scatter = ax.scatter( + _r, _l, s=100, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7 + ) # Plot horizontal line ax.axhline(y=1, color="k", linestyle="--") @@ -97,30 +103,67 @@ cax = divider.append_axes("right", size="5%", pad=0.05) cax.tick_params(labelsize=14) fig.colorbar(scatter, cax=cax) -ax.text(1e7, 0.013, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=12) +ax.text( + 1e7, + 0.013, + r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", + fontsize=12, +) # Legend... -blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None', - markersize=8, label=r'$n_x = 960$') -red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None', - markersize=8, label=r'$n_x = 1920$') -purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 3840$') -diamond = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 7680$') - -ax.legend(handles=[blue_star, red_square, purple_triangle, diamond], - loc="upper center", - bbox_to_anchor=(0.5,1.1), - borderaxespad=0., - ncol=len(markers), - handletextpad=0.1, - fontsize=14) +blue_star = mlines.Line2D( + [], + [], + color="red", + marker="o", + linestyle="None", + markersize=8, + label=r"$n_x = 960$", +) +red_square = mlines.Line2D( + [], + [], + color="red", + marker="s", + linestyle="None", + markersize=8, + label=r"$n_x = 1920$", +) +purple_triangle = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 3840$", +) +diamond = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 7680$", +) + +ax.legend( + handles=[blue_star, red_square, purple_triangle, diamond], + loc="upper center", + bbox_to_anchor=(0.5, 1.1), + borderaxespad=0.0, + ncol=len(markers), + handletextpad=0.1, + fontsize=14, +) SAVE = False if SAVE: - path_save = "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + \ - "vertical_length_scale.png" + path_save = ( + "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + + "vertical_length_scale.png" + ) fig.savefig(path_save, format="png", bbox_inches="tight") plt.show() diff --git a/scripts/ns2d.strat/plot_vertical_scale_godoy.py b/scripts/ns2d.strat/plot_vertical_scale_godoy.py index d8995af8..5f07e0f4 100644 --- a/scripts/ns2d.strat/plot_vertical_scale_godoy.py +++ b/scripts/ns2d.strat/plot_vertical_scale_godoy.py @@ -20,12 +20,14 @@ from fluiddyn.output.rcparams import set_rcparams path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim" -directories = ["sim960_no_shear_modes", - "sim960_no_shear_modes_transitory", - "sim1920_no_shear_modes", - "sim1920_modif_res_no_shear_modes", - "sim3840_modif_res_no_shear_modes", - "sim7680_modif_res_no_shear_modes"] +directories = [ + "sim960_no_shear_modes", + "sim960_no_shear_modes_transitory", + "sim1920_no_shear_modes", + "sim1920_modif_res_no_shear_modes", + "sim3840_modif_res_no_shear_modes", + "sim7680_modif_res_no_shear_modes", +] paths_simulations = [] @@ -68,12 +70,12 @@ for path_file in path_phys_files[-10:-1]: with h5py.File(path_file, "r") as f: ux = f["state_phys"]["ux"][...] - ux_rms.append(np.sqrt(np.mean(ux**2))) + ux_rms.append(np.sqrt(np.mean(ux ** 2))) # Load parameters params = load_params_simul(path) - compensate_godoy = l_x / (Re_8**(1/8)) + compensate_godoy = l_x / (Re_8 ** (1 / 8)) lzs_godoy.append(l_z / compensate_godoy) # Parameters figures @@ -87,8 +89,12 @@ ax.set_ylim([0, 1.2e1]) # Plot... -for _f, _r, _a, _d, _m, _l in zip(froudes, reynoldsb, anisotropies, dissipations, markers, lzs_godoy): - scatter = ax.scatter(_r, _l, s=100, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7) +for _f, _r, _a, _d, _m, _l in zip( + froudes, reynoldsb, anisotropies, dissipations, markers, lzs_godoy +): + scatter = ax.scatter( + _r, _l, s=100, c=_d, vmin=0, vmax=1, marker=_m, alpha=0.7 + ) # Plot horizontal line ax.axhline(y=1, color="k", linestyle="--") @@ -98,30 +104,67 @@ cax = divider.append_axes("right", size="5%", pad=0.05) cax.tick_params(labelsize=14) fig.colorbar(scatter, cax=cax) -ax.text(1e7, 0.013, r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", fontsize=12) +ax.text( + 1e7, + 0.013, + r"$\log_{10} \left(\frac{k_{x, 1/2}}{k_{x, f}}\right)$", + fontsize=12, +) # Legend... -blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None', - markersize=8, label=r'$n_x = 960$') -red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None', - markersize=8, label=r'$n_x = 1920$') -purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 3840$') -diamond = mlines.Line2D([], [], color='red', marker='^', linestyle='None', - markersize=8, label=r'$n_x = 7680$') - -ax.legend(handles=[blue_star, red_square, purple_triangle, diamond], - loc="upper center", - bbox_to_anchor=(0.5,1.1), - borderaxespad=0., - ncol=len(markers), - handletextpad=0.1, - fontsize=14) +blue_star = mlines.Line2D( + [], + [], + color="red", + marker="o", + linestyle="None", + markersize=8, + label=r"$n_x = 960$", +) +red_square = mlines.Line2D( + [], + [], + color="red", + marker="s", + linestyle="None", + markersize=8, + label=r"$n_x = 1920$", +) +purple_triangle = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 3840$", +) +diamond = mlines.Line2D( + [], + [], + color="red", + marker="^", + linestyle="None", + markersize=8, + label=r"$n_x = 7680$", +) + +ax.legend( + handles=[blue_star, red_square, purple_triangle, diamond], + loc="upper center", + bbox_to_anchor=(0.5, 1.1), + borderaxespad=0.0, + ncol=len(markers), + handletextpad=0.1, + fontsize=14, +) SAVE = True if SAVE: - path_save = "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + \ - "vertical_length_scale_godoy.png" + path_save = ( + "/fsnet/project/meige/2015/15DELDUCA/notebooks/figures/" + + "vertical_length_scale_godoy.png" + ) fig.savefig(path_save, format="png", bbox_inches="tight") plt.show() diff --git a/scripts/ns2d.strat/simul_from_state.py b/scripts/ns2d.strat/simul_from_state.py index 1373bc6c..cf858213 100644 --- a/scripts/ns2d.strat/simul_from_state.py +++ b/scripts/ns2d.strat/simul_from_state.py @@ -57,7 +57,7 @@ params.output.periods_save.spatial_means = 0.5 params.output.periods_save.spectra = 0.5 params.output.periods_save.spect_energy_budg = 0.5 -params.output.periods_save.spatio_temporal_spectra = 1. +params.output.periods_save.spatio_temporal_spectra = 1.0 params.output.spatio_temporal_spectra.size_max_file = 100 params.output.spatio_temporal_spectra.spatial_decimate = 4 diff --git a/scripts/ns2d.strat/simul_init_linear_mode.py b/scripts/ns2d.strat/simul_init_linear_mode.py index f8013af2..db7f0e3c 100644 --- a/scripts/ns2d.strat/simul_init_linear_mode.py +++ b/scripts/ns2d.strat/simul_init_linear_mode.py @@ -22,18 +22,18 @@ params = Simul.create_default_params() -params.N = 50. +params.N = 50.0 params.oper.nx = nx = 128 params.oper.ny = ny = nx // 4 # Parameters time stepping -params.time_stepping.USE_CFL = True -params.time_stepping.t_end = 2. +params.time_stepping.USE_CFL = True +params.time_stepping.t_end = 2.0 # Field initialization in the script params.init_fields.type = "linear_mode" params.init_fields.linear_mode.eigenmode = "ap_fft" -params.init_fields.linear_mode.i_mode = (2,1) +params.init_fields.linear_mode.i_mode = (2, 1) params.init_fields.linear_mode.delta_k_adim = 1 # Parameters output diff --git a/scripts/ns2d.strat/submit_legi.py b/scripts/ns2d.strat/submit_legi.py index 86adc5ad..3d23bc2c 100644 --- a/scripts/ns2d.strat/submit_legi.py +++ b/scripts/ns2d.strat/submit_legi.py @@ -10,11 +10,12 @@ cluster = Calcul() cluster.commands_setting_env.append( - 'export FLUIDSIM_PATH="/fsnet/project/meige/2015/15DELDUCA/DataSim"') + 'export FLUIDSIM_PATH="/fsnet/project/meige/2015/15DELDUCA/DataSim"' +) name_run_root = f"find_coeff_nu8_gamma={args.gamma}" -walltime = '24:00:00' +walltime = "24:00:00" nb_proc = 8 command_to_submit = f"python coeff_diss.py {args.gamma}" @@ -26,4 +27,7 @@ walltime=walltime, nb_mpi_processes=nb_proc, omp_num_threads=1, - idempotent=True, delay_signal_walltime=300, ask=False) + idempotent=True, + delay_signal_walltime=300, + ask=False, +) diff --git a/scripts/ns2d.strat/submit_simul_from_state.py b/scripts/ns2d.strat/submit_simul_from_state.py index baf460c1..557f04e9 100644 --- a/scripts/ns2d.strat/submit_simul_from_state.py +++ b/scripts/ns2d.strat/submit_simul_from_state.py @@ -12,17 +12,20 @@ cluster = Calcul() cluster.commands_setting_env.append( - 'export FLUIDSIM_PATH="/fsnet/project/meige/2015/15DELDUCA/DataSim"') + 'export FLUIDSIM_PATH="/fsnet/project/meige/2015/15DELDUCA/DataSim"' +) name_run_root = "sim480_gamma={}_NO_SHEAR_MODES={}".format( - args.gamma, bool(args.NO_SHEAR_MODES)) + args.gamma, bool(args.NO_SHEAR_MODES) +) -walltime = '24:00:00' +walltime = "24:00:00" # nb_proc = cluster.nb_cores_per_node nb_proc = 8 command_to_submit = "python simul_from_state.py {} {}".format( - args.gamma, args.NO_SHEAR_MODES) + args.gamma, args.NO_SHEAR_MODES +) cluster.submit_command( command_to_submit, @@ -31,4 +34,7 @@ walltime=walltime, nb_mpi_processes=nb_proc, omp_num_threads=1, - idempotent=True, delay_signal_walltime=300, ask=False) + idempotent=True, + delay_signal_walltime=300, + ask=False, +) diff --git a/scripts/ns2d.strat/tests/check_freq_spectra_with_linear_mode.py b/scripts/ns2d.strat/tests/check_freq_spectra_with_linear_mode.py index e30c4c4a..d6cecea1 100644 --- a/scripts/ns2d.strat/tests/check_freq_spectra_with_linear_mode.py +++ b/scripts/ns2d.strat/tests/check_freq_spectra_with_linear_mode.py @@ -27,6 +27,7 @@ from glob import glob from fluidsim.solvers.ns2d.strat.solver import Simul + def _create_object_params(): params = Simul.create_default_params() try: @@ -40,10 +41,10 @@ def _create_object_params(): # Forcing parameters params.forcing.enable = False - params.forcing.type = 'tcrandom_anisotropic' + params.forcing.type = "tcrandom_anisotropic" try: - params.forcing.tcrandom_anisotropic.angle = '45.0°' + params.forcing.tcrandom_anisotropic.angle = "45.0°" except AttributeError: pass @@ -52,6 +53,7 @@ def _create_object_params(): # Compute \omega_l from math import radians + if "°" in params.forcing.tcrandom_anisotropic.angle: angle = params.forcing.tcrandom_anisotropic.angle.split("°") angle = float(angle[0]) @@ -60,21 +62,21 @@ def _create_object_params(): omega_l = params.N * np.sin(radians(angle)) params.forcing.tcrandom.time_correlation = 2 * pi / omega_l - params.forcing.key_forced = 'ap_fft' + params.forcing.key_forced = "ap_fft" # Time stepping parameters params.time_stepping.USE_CFL = True params.time_stepping.USE_T_END = True - params.time_stepping.t_end = 2. + params.time_stepping.t_end = 2.0 # Output parameters params.output.HAS_TO_SAVE = False - params.output.sub_directory = 'tests' + params.output.sub_directory = "tests" return params -if __name__ == '__main__': +if __name__ == "__main__": SAVE = True format = ".pdf" ### SHORT SIMULATION ### @@ -87,22 +89,22 @@ def _create_object_params(): params.oper.Ly = params.oper.Lx * (ny / nx) params.oper.NO_SHEAR_MODES = False - params.nu_8 = 0. - params.N = 50. + params.nu_8 = 0.0 + params.N = 50.0 params.time_stepping.USE_CFL = False params.time_stepping.deltat0 = 0.005 - params.time_stepping.t_end = 10. + params.time_stepping.t_end = 10.0 params.time_stepping.cfl_coef_group = None params.output.HAS_TO_SAVE = True - params.output.periods_print.print_stdout = 1. + params.output.periods_print.print_stdout = 1.0 params.output.periods_save.phys_fields = 2e-1 params.output.periods_save.spatial_means = 0.0005 params.output.periods_save.frequency_spectra = 1 - params.output.frequency_spectra.time_start = 0. + params.output.frequency_spectra.time_start = 0.0 params.output.frequency_spectra.spatial_decimate = 1 params.output.frequency_spectra.size_max_file = 10 params.output.frequency_spectra.time_decimate = 4 @@ -117,14 +119,19 @@ def _create_object_params(): sim.time_stepping.start() from fluiddyn.util import mpi + if mpi.rank == 0: kx_s = sim.oper.KX[params.init_fields.linear_mode.i_mode] kz_s = sim.oper.KY[params.init_fields.linear_mode.i_mode] from math import pi - omega_n = params.N * np.sin(np.arctan( - sim.oper.kx[params.init_fields.linear_mode.i_mode[0]]/ \ - sim.oper.ky[params.init_fields.linear_mode.i_mode[1]])) + + omega_n = params.N * np.sin( + np.arctan( + sim.oper.kx[params.init_fields.linear_mode.i_mode[0]] + / sim.oper.ky[params.init_fields.linear_mode.i_mode[1]] + ) + ) omega_n = omega_n / (2 * pi) @@ -132,14 +139,16 @@ def _create_object_params(): sim.output.frequency_spectra.compute_frequency_spectra() ### LOAD DATA AND PLOT ### - path_file = glob(os.path.join( - sim.output.path_run, "temporal_data", "temp_*"))[0] + path_file = glob( + os.path.join(sim.output.path_run, "temporal_data", "temp_*") + )[0] with h5py.File(path_file, "r") as f: omegas = f["omegas"][...] freq_spectrum = f["freq_spectrum"][...] import matplotlib.pyplot as plt + fig, ax = plt.subplots() ax.set_xlabel(r"$\omega / \omega_{th}$", fontsize=16) ax.set_ylabel(r"F($\omega$)", fontsize=16) @@ -148,13 +157,16 @@ def _create_object_params(): # for i in range(0,10): # ax.loglog(omegas/omega_n, freq_spectrum[0, :, 1, i]) - ax.semilogy(omegas/omega_n, freq_spectrum[0, :, 1, 4]) + ax.semilogy(omegas / omega_n, freq_spectrum[0, :, 1, 4]) ax.axvline(x=omega_n / omega_n, label=r"$\omega_{th}$", c="k") # Set text - ax.text(5e-2, 1e3, - r"$\omega_{th} = N \sin(arctan \left( \frac{k_x}{k_z} \right))$", - fontsize=16) + ax.text( + 5e-2, + 1e3, + r"$\omega_{th} = N \sin(arctan \left( \frac{k_x}{k_z} \right))$", + fontsize=16, + ) # If SAVE: if SAVE: @@ -162,7 +174,10 @@ def _create_object_params(): path_save = path_root_save + f"/test_frequency_spectra_seq{format}" if mpi.nb_proc > 1: - path_save = path_root_save + f"/test_frequency_spectra_mpi_{mpi.nb_proc}{format}" + path_save = ( + path_root_save + + f"/test_frequency_spectra_mpi_{mpi.nb_proc}{format}" + ) fig.savefig(path_save, format="pdf") ax.legend() diff --git a/scripts/plot_results/dispersion_relationship_SW1l.py b/scripts/plot_results/dispersion_relationship_SW1l.py index d5665c80..cc3c587e 100644 --- a/scripts/plot_results/dispersion_relationship_SW1l.py +++ b/scripts/plot_results/dispersion_relationship_SW1l.py @@ -6,51 +6,41 @@ import matplotlib.pyplot as plt -f = 1. -c= 1. +f = 1.0 +c = 1.0 -L_D = c/f +L_D = c / f -Lx = 30. +Lx = 30.0 nx = 64 -nkx = nx/2. +1 +nkx = nx / 2.0 + 1 -kx = 2*np.pi/Lx* np.arange(nkx) +kx = 2 * np.pi / Lx * np.arange(nkx) -omega_theo = np.sqrt( f**2 + (c*kx)**2 ) +omega_theo = np.sqrt(f ** 2 + (c * kx) ** 2) # numerical results -k_num = np.array([0.63, 1.88, 2.51, 4.19 -]) - -frequency = np.array([0.1875, 0.34, 0.43, 0.7 -]) - - +k_num = np.array([0.63, 1.88, 2.51, 4.19]) +frequency = np.array([0.1875, 0.34, 0.43, 0.7]) ikx = 14 - - - fig = plt.figure(2) fig.clf() plt.hold(True) -plt.plot(kx*L_D, omega_theo/(2*np.pi), 'k', linewidth=2 ) +plt.plot(kx * L_D, omega_theo / (2 * np.pi), "k", linewidth=2) -plt.plot(kx*L_D, kx*L_D/(2*np.pi), 'k--', linewidth=0.5 ) +plt.plot(kx * L_D, kx * L_D / (2 * np.pi), "k--", linewidth=0.5) -plt.plot(kx[ikx]*L_D, omega_theo[ikx]/(2*np.pi), 'rx', linewidth=2 ) +plt.plot(kx[ikx] * L_D, omega_theo[ikx] / (2 * np.pi), "rx", linewidth=2) -plt.plot(k_num*L_D, frequency, 'bo', linewidth=2 ) +plt.plot(k_num * L_D, frequency, "bo", linewidth=2) -plt.xlabel('k L_D') -plt.ylabel('frequency') +plt.xlabel("k L_D") +plt.ylabel("frequency") plt.show() - - diff --git a/scripts/plot_results/fig_spectra_forcingq.py b/scripts/plot_results/fig_spectra_forcingq.py index 7a4ab41c..15c9862c 100644 --- a/scripts/plot_results/fig_spectra_forcingq.py +++ b/scripts/plot_results/fig_spectra_forcingq.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -#coding=utf8 +# coding=utf8 import numpy as np import matplotlib.pyplot as plt @@ -11,109 +11,83 @@ num_fig = 1000 SAVE_FIG = 1 -name_file = 'spectra_forcingq.eps' +name_file = "spectra_forcingq.eps" create_fig = CreateFigArticles( - short_name_article='SW1l', - SAVE_FIG=SAVE_FIG, - FOR_BEAMER=False, - fontsize=19 - ) - - + short_name_article="SW1l", SAVE_FIG=SAVE_FIG, FOR_BEAMER=False, fontsize=19 +) name_dir_results = ( -create_fig.path_base_dir+'/Results_for_article_SW1l' -'/Approach_runs_2048x2048' -'/SE2D_SW1lexlin_forcing_L=50.x50._2048x2048_c2=400_f=0_2013-05-29_23-54-57' + create_fig.path_base_dir + "/Results_for_article_SW1l" + "/Approach_runs_2048x2048" + "/SE2D_SW1lexlin_forcing_L=50.x50._2048x2048_c2=400_f=0_2013-05-29_23-54-57" ) sim = solveq2d.create_sim_plot_from_dir(name_dir_results) tmin = 30 dict_results = sim.output.spectra.load2D_mean(tmin=tmin) -kh = dict_results['kh'] +kh = dict_results["kh"] -EK = dict_results['spectrum2D_EK'] -EA = dict_results['spectrum2D_EA'] -EKr = dict_results['spectrum2D_EKr'] +EK = dict_results["spectrum2D_EK"] +EA = dict_results["spectrum2D_EA"] +EKr = dict_results["spectrum2D_EKr"] E_tot = EK + EA EKd = EK - EKr -Edlin = dict_results['spectrum2D_Edlin'] +Edlin = dict_results["spectrum2D_Edlin"] fig, ax1 = create_fig.figure_axe(name_file=name_file) -ax1.set_xscale('log') -ax1.set_yscale('log') - -coef_compensate = 5./3 -coef_norm = kh**coef_compensate - -l_Etot = ax1.plot(kh, E_tot*coef_norm, 'k', linewidth=4) -l_EK = ax1.plot(kh, EK*coef_norm, 'r', linewidth=2) -l_EA = ax1.plot(kh, EA*coef_norm, 'b', linewidth=2) -ax1.plot(kh, EKr*coef_norm, 'r--', linewidth=2) -ax1.plot(kh, EKd*coef_norm, 'r:', linewidth=2) - -ax1.plot(kh, -EK*coef_norm, 'm', linewidth=2) -ax1.plot(kh, -EKd*coef_norm, 'm:', linewidth=2) - -ax1.plot(kh, Edlin*coef_norm, 'y:', linewidth=2) - - -cond = np.logical_and(kh > 1 , kh < 20) -ax1.plot(kh[cond], 1e1*kh[cond]**(-3.)*coef_norm[cond], 'k--', linewidth=1) -plt.figtext(0.6, 0.78, '$k^{-3}$', fontsize=20) - -cond = np.logical_and(kh > 0.3 , kh < 10) -ax1.plot(kh[cond], 4e-2*kh[cond]**(-2.)*coef_norm[cond], 'k:', linewidth=1) -plt.figtext(0.25, 0.55, '$k^{-2}$', fontsize=20) - -cond = np.logical_and(kh > 0.3 , kh < 15) -ax1.plot(kh[cond], 1e-2*kh[cond]**(-5./3)*coef_norm[cond], - 'k-.', linewidth=1) -plt.figtext(0.5, 0.45, '$k^{-5/3}$', fontsize=20) - - +ax1.set_xscale("log") +ax1.set_yscale("log") +coef_compensate = 5.0 / 3 +coef_norm = kh ** coef_compensate +l_Etot = ax1.plot(kh, E_tot * coef_norm, "k", linewidth=4) +l_EK = ax1.plot(kh, EK * coef_norm, "r", linewidth=2) +l_EA = ax1.plot(kh, EA * coef_norm, "b", linewidth=2) +ax1.plot(kh, EKr * coef_norm, "r--", linewidth=2) +ax1.plot(kh, EKd * coef_norm, "r:", linewidth=2) +ax1.plot(kh, -EK * coef_norm, "m", linewidth=2) +ax1.plot(kh, -EKd * coef_norm, "m:", linewidth=2) +ax1.plot(kh, Edlin * coef_norm, "y:", linewidth=2) +cond = np.logical_and(kh > 1, kh < 20) +ax1.plot(kh[cond], 1e1 * kh[cond] ** (-3.0) * coef_norm[cond], "k--", linewidth=1) +plt.figtext(0.6, 0.78, "$k^{-3}$", fontsize=20) +cond = np.logical_and(kh > 0.3, kh < 10) +ax1.plot(kh[cond], 4e-2 * kh[cond] ** (-2.0) * coef_norm[cond], "k:", linewidth=1) +plt.figtext(0.25, 0.55, "$k^{-2}$", fontsize=20) -ax1.set_xlabel('$k_h$') -ax1.set_ylabel('2D spectra') +cond = np.logical_and(kh > 0.3, kh < 15) +ax1.plot( + kh[cond], 1e-2 * kh[cond] ** (-5.0 / 3) * coef_norm[cond], "k-.", linewidth=1 +) +plt.figtext(0.5, 0.45, "$k^{-5/3}$", fontsize=20) +ax1.set_xlabel("$k_h$") +ax1.set_ylabel("2D spectra") -plt.rc('legend', numpoints=1) +plt.rc("legend", numpoints=1) leg1 = plt.figlegend( - [l_Etot[0], l_EK[0], l_EA[0]], - ['$E$', '$E_K$', '$E_A$'], - loc=(0.78, 0.7), - labelspacing = 0.2 + [l_Etot[0], l_EK[0], l_EA[0]], + ["$E$", "$E_K$", "$E_A$"], + loc=(0.78, 0.7), + labelspacing=0.2, ) - -ax1.set_xlim([0.1,150]) -ax1.set_ylim([1e-5,2e1]) +ax1.set_xlim([0.1, 150]) +ax1.set_ylim([1e-5, 2e1]) create_fig.save_fig() plt.show() - - - - - - - - - - - diff --git a/scripts/plot_results/fig_spectra_forcingq_diff_c.py b/scripts/plot_results/fig_spectra_forcingq_diff_c.py index f33f400f..89d772e3 100644 --- a/scripts/plot_results/fig_spectra_forcingq_diff_c.py +++ b/scripts/plot_results/fig_spectra_forcingq_diff_c.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -#coding=utf8 +# coding=utf8 import numpy as np import matplotlib.pyplot as plt @@ -11,20 +11,18 @@ num_fig = 1000 SAVE_FIG = 1 -name_file = 'spectra_forcingq_diff_c.eps' +name_file = "spectra_forcingq_diff_c.eps" create_fig = CreateFigArticles( - short_name_article='SW1l', - SAVE_FIG=SAVE_FIG, - FOR_BEAMER=False, - fontsize=19) + short_name_article="SW1l", SAVE_FIG=SAVE_FIG, FOR_BEAMER=False, fontsize=19 +) dir_base = ( - create_fig.path_base_dir+'/Results_for_article_SW1l' - '/Approach_runs_2048x2048') + create_fig.path_base_dir + "/Results_for_article_SW1l" + "/Approach_runs_2048x2048" +) set_of_dir_results = solveq2d.SetOfDirResults(dir_base=dir_base) -dirs = set_of_dir_results.dirs_from_values(solver='SW1lexlin', - FORCING=True) +dirs = set_of_dir_results.dirs_from_values(solver="SW1lexlin", FORCING=True) print(dirs) @@ -35,76 +33,57 @@ def sprectra_from_namedir(name_dir_results): path_dir_results = set_of_dir_results.path_dirs[name_dir_results] sim = solveq2d.create_sim_plot_from_dir(path_dir_results) dict_results = sim.output.spectra.load2D_mean(tmin=tmin) - kh = dict_results['kh'] - EK = dict_results['spectrum2D_EK'] - EA = dict_results['spectrum2D_EA'] - EKr = dict_results['spectrum2D_EKr'] + kh = dict_results["kh"] + EK = dict_results["spectrum2D_EK"] + EA = dict_results["spectrum2D_EA"] + EKr = dict_results["spectrum2D_EKr"] EKd = EK - EKr return kh, EKr, EKd - kh, EKr0, EKd0 = sprectra_from_namedir(dirs[0]) kh, EKr1, EKd1 = sprectra_from_namedir(dirs[1]) kh, EKr2, EKd2 = sprectra_from_namedir(dirs[2]) kh, EKr3, EKd3 = sprectra_from_namedir(dirs[3]) - - - - fig, ax1 = create_fig.figure_axe(name_file=name_file) -ax1.set_xscale('log') -ax1.set_yscale('log') - -coef_compensate = 5./3 -coef_norm = kh**coef_compensate - -ax1.plot(kh, EKr0*coef_norm, 'k--', linewidth=2) -ax1.plot(kh, EKd0*coef_norm, 'k:', linewidth=2) - -ax1.plot(kh, EKr1*coef_norm, 'r--', linewidth=2) -ax1.plot(kh, EKd1*coef_norm, 'r:', linewidth=2) - -ax1.plot(kh, EKr2*coef_norm, 'b--', linewidth=2) -ax1.plot(kh, EKd2*coef_norm, 'b:', linewidth=2) - -ax1.plot(kh, EKr3*coef_norm, 'c--', linewidth=2) -ax1.plot(kh, EKd3*coef_norm, 'c:', linewidth=2) +ax1.set_xscale("log") +ax1.set_yscale("log") +coef_compensate = 5.0 / 3 +coef_norm = kh ** coef_compensate +ax1.plot(kh, EKr0 * coef_norm, "k--", linewidth=2) +ax1.plot(kh, EKd0 * coef_norm, "k:", linewidth=2) +ax1.plot(kh, EKr1 * coef_norm, "r--", linewidth=2) +ax1.plot(kh, EKd1 * coef_norm, "r:", linewidth=2) +ax1.plot(kh, EKr2 * coef_norm, "b--", linewidth=2) +ax1.plot(kh, EKd2 * coef_norm, "b:", linewidth=2) -cond = np.logical_and(kh > 1 , kh < 20) -ax1.plot(kh[cond], 1e1*kh[cond]**(-3.)*coef_norm[cond], - 'k--', linewidth=1) -plt.figtext(0.6, 0.78, '$k^{-3}$', fontsize=20) +ax1.plot(kh, EKr3 * coef_norm, "c--", linewidth=2) +ax1.plot(kh, EKd3 * coef_norm, "c:", linewidth=2) -cond = np.logical_and(kh > 0.3 , kh < 10) -ax1.plot(kh[cond], 4e-2*kh[cond]**(-2.)*coef_norm[cond], - 'k:', linewidth=1) -plt.figtext(0.25, 0.55, '$k^{-2}$', fontsize=20) -cond = np.logical_and(kh > 0.3 , kh < 15) -ax1.plot(kh[cond], 1e-2*kh[cond]**(-5./3)*coef_norm[cond], - 'k-.', linewidth=1) -plt.figtext(0.5, 0.45, '$k^{-5/3}$', fontsize=20) +cond = np.logical_and(kh > 1, kh < 20) +ax1.plot(kh[cond], 1e1 * kh[cond] ** (-3.0) * coef_norm[cond], "k--", linewidth=1) +plt.figtext(0.6, 0.78, "$k^{-3}$", fontsize=20) +cond = np.logical_and(kh > 0.3, kh < 10) +ax1.plot(kh[cond], 4e-2 * kh[cond] ** (-2.0) * coef_norm[cond], "k:", linewidth=1) +plt.figtext(0.25, 0.55, "$k^{-2}$", fontsize=20) +cond = np.logical_and(kh > 0.3, kh < 15) +ax1.plot( + kh[cond], 1e-2 * kh[cond] ** (-5.0 / 3) * coef_norm[cond], "k-.", linewidth=1 +) +plt.figtext(0.5, 0.45, "$k^{-5/3}$", fontsize=20) - - - - - - -ax1.set_xlabel('$k_h$') -ax1.set_ylabel('2D spectra') - - +ax1.set_xlabel("$k_h$") +ax1.set_ylabel("2D spectra") # plt.rc('legend', numpoints=1) @@ -116,9 +95,8 @@ def sprectra_from_namedir(name_dir_results): # ) - -ax1.set_xlim([0.1,150]) -ax1.set_ylim([1e-5,2e1]) +ax1.set_xlim([0.1, 150]) +ax1.set_ylim([1e-5, 2e1]) create_fig.save_fig() diff --git a/scripts/plot_results/plot_interm.py b/scripts/plot_results/plot_interm.py index e9167785..73d72768 100644 --- a/scripts/plot_results/plot_interm.py +++ b/scripts/plot_results/plot_interm.py @@ -1,5 +1,3 @@ - - import matplotlib.pylab as plt import glob import numpy as np @@ -13,21 +11,19 @@ SAVE_FIG = False create_fig = CreateFigArticles( - short_name_article='SW1l', - SAVE_FIG=SAVE_FIG, - FOR_BEAMER=False, - fontsize=19) + short_name_article="SW1l", SAVE_FIG=SAVE_FIG, FOR_BEAMER=False, fontsize=19 +) -dir_base = create_fig.path_base_dir+'/Results_SW1lw' +dir_base = create_fig.path_base_dir + "/Results_SW1lw" c = 40 -resol = 240*2**5 +resol = 240 * 2 ** 5 str_resol = repr(resol) str_to_find_path = ( - dir_base+'/Pure_standing_waves_'+ - str_resol+'*/SE2D*c='+repr(c))+'_*' + dir_base + "/Pure_standing_waves_" + str_resol + "*/SE2D*c=" + repr(c) +) + "_*" print(str_to_find_path) paths_dir = glob.glob(str_to_find_path) @@ -38,33 +34,30 @@ sim = solveq2d.create_sim_plot_from_dir(paths_dir[0]) tmin = sim.output.spatial_means.first_saved_time() -tstatio = tmin + 4. - +tstatio = tmin + 4.0 tmax = 1000 +key_var = "uy" - -key_var = 'uy' - -(pdf_timemean, values_inc_timemean, nb_rx_to_plot - ) = sim.output.increments.load_pdf_from_file( - tmin=tmin, tmax=tmax, key_var=key_var) - - - - +( + pdf_timemean, + values_inc_timemean, + nb_rx_to_plot, +) = sim.output.increments.load_pdf_from_file( + tmin=tmin, tmax=tmax, key_var=key_var +) -deltax = sim.param.Lx/sim.param.nx -rxs = np.array(sim.output.increments.rxs, dtype=np.float64)*deltax +deltax = sim.param.Lx / sim.param.nx +rxs = np.array(sim.output.increments.rxs, dtype=np.float64) * deltax # if 7680 -rmin = 8*deltax -rmax = 40*deltax +rmin = 8 * deltax +rmax = 40 * deltax # rmin = 0.6 # rmax = 3 @@ -83,10 +76,7 @@ # rmax = 3e-1 - - -condr = np.logical_and(rxs>rmin, rxs rmin, rxs < rmax) def expo_from_order(order, PLOT=False, PLOT_PDF=False): @@ -94,9 +84,9 @@ def expo_from_order(order, PLOT=False, PLOT_PDF=False): M_order = np.empty(rxs.shape) for irx in xrange(rxs.size): deltainc = values_inc_timemean[irx, 1] - values_inc_timemean[irx, 0] - M_order[irx] = deltainc*np.sum( - pdf_timemean[irx] - *abs(values_inc_timemean[irx])**order) + M_order[irx] = deltainc * np.sum( + pdf_timemean[irx] * abs(values_inc_timemean[irx]) ** order + ) # M_order[irx] = np.abs(deltainc*np.sum( # pdf_timemean[irx] @@ -106,90 +96,93 @@ def expo_from_order(order, PLOT=False, PLOT_PDF=False): pol = np.polyfit(np.log(rxs[condr]), np.log(M_order[condr]), 1) expo = pol[0] - print(f'order = {order:.2f} ; expo = {expo:.2f}') - M_lin = np.exp((pol[1] + np.log(rxs[condr])*pol[0])) - + print(f"order = {order:.2f} ; expo = {expo:.2f}") + M_lin = np.exp((pol[1] + np.log(rxs[condr]) * pol[0])) if PLOT: fig, ax1 = sim.output.figure_axe() title = ( - 'struct. function, solver '+sim.output.name_solver+ - f', nh = {sim.param.nx:5d}'+ - ', c = {:.4g}, f = {:.4g}'.format(np.sqrt(sim.param.c2), - sim.param.f) + "struct. function, solver " + + sim.output.name_solver + + f", nh = {sim.param.nx:5d}" + + ", c = {:.4g}, f = {:.4g}".format( + np.sqrt(sim.param.c2), sim.param.f ) + ) - ax1.set_xlabel('$r$') - ax1.set_ylabel( - r'$\langle |\delta v|^{'+f'{order:.2f}'+'}\\rangle/r$') + ax1.set_xlabel("$r$") + ax1.set_ylabel(r"$\langle |\delta v|^{" + f"{order:.2f}" + "}\\rangle/r$") ax1.set_title(title) ax1.hold(True) - ax1.set_xscale('log') - ax1.set_yscale('log') + ax1.set_xscale("log") + ax1.set_yscale("log") - norm = rxs**(1) + norm = rxs ** (1) - ax1.plot(rxs, M_order/norm, - 'x-', linewidth=1) + ax1.plot(rxs, M_order / norm, "x-", linewidth=1) - ax1.plot(rxs[condr], M_order[condr]/norm[condr], 'x-r', linewidth=1) + ax1.plot(rxs[condr], M_order[condr] / norm[condr], "x-r", linewidth=1) - ax1.plot(rxs[condr], M_lin/norm[condr], 'y', linewidth=2) + ax1.plot(rxs[condr], M_lin / norm[condr], "y", linewidth=2) - l_1 = ax1.plot(rxs, rxs/norm, 'k', linewidth=1) - l_K41 = ax1.plot(rxs, rxs**(order/3)/norm, 'k--', linewidth=1) + l_1 = ax1.plot(rxs, rxs / norm, "k", linewidth=1) + l_K41 = ax1.plot(rxs, rxs ** (order / 3) / norm, "k--", linewidth=1) - cond = rxs < 4*deltax - temp = rxs**(order)/norm - l_smooth= ax1.plot( - rxs[cond],temp[cond]/temp[0]*M_order[0]/norm[0], - 'k:', linewidth=1) + cond = rxs < 4 * deltax + temp = rxs ** (order) / norm + l_smooth = ax1.plot( + rxs[cond], + temp[cond] / temp[0] * M_order[0] / norm[0], + "k:", + linewidth=1, + ) - plt.text(1, 6*10**0,f'order = {order:.2f}',fontsize=16) - plt.text(0.11, 4*10**-1,f'expo = {expo:.3f}',fontsize=16) + plt.text(1, 6 * 10 ** 0, f"order = {order:.2f}", fontsize=16) + plt.text(0.11, 4 * 10 ** -1, f"expo = {expo:.3f}", fontsize=16) - - - plt.rc('legend', numpoints=1) + plt.rc("legend", numpoints=1) leg1 = ax1.legend( [l_smooth[0], l_K41[0], l_1[0]], - ['smooth $r^q$', 'K41 $r^{q/3}$', - 'shocks $r^1$'], + ["smooth $r^q$", "K41 $r^{q/3}$", "shocks $r^1$"], loc=0 # labelspacing = 0.2 -) - - - + ) if PLOT_PDF: fig, ax1 = sim.output.figure_axe() - title = ('pdf increments, solver '+sim.output.name_solver+ -f', nh = {sim.param.nx:5d}'+ -f', c2 = {sim.param.c2:.4g}, f = {sim.param.f:.4g}' - ) + title = ( + "pdf increments, solver " + + sim.output.name_solver + + f", nh = {sim.param.nx:5d}" + + f", c2 = {sim.param.c2:.4g}, f = {sim.param.f:.4g}" + ) ax1.set_title(title) ax1.hold(True) - ax1.set_xscale('linear') - ax1.set_yscale('linear') + ax1.set_xscale("linear") + ax1.set_yscale("linear") ax1.set_xlabel(key_var) - ax1.set_ylabel(r'PDF x $\delta v^'+repr(order)+'$') + ax1.set_ylabel(r"PDF x $\delta v^" + repr(order) + "$") - colors = ['k', 'y', 'r', 'b', 'g', 'm', 'c'] + colors = ["k", "y", "r", "b", "g", "m", "c"] irx_to_plot = [10, 50, 100] for irxp, irx in enumerate(irx_to_plot): val_inc = values_inc_timemean[irx] - ax1.plot(val_inc, pdf_timemean[irx]*abs(val_inc)**order, - colors[irxp]+'x-', linewidth=1) + ax1.plot( + val_inc, + pdf_timemean[irx] * abs(val_inc) ** order, + colors[irxp] + "x-", + linewidth=1, + ) return expo + PLOT = False PLOT = True if PLOT: @@ -198,31 +191,33 @@ def expo_from_order(order, PLOT=False, PLOT_PDF=False): delta_order = 0.05 -orders = np.arange(0., 6, delta_order) +orders = np.arange(0.0, 6, delta_order) expos = np.empty(orders.shape) for iorder, order in enumerate(orders): expos[iorder] = expo_from_order(order, PLOT=PLOT, PLOT_PDF=False) -expos_K41 = orders/3 +expos_K41 = orders / 3 fig, ax1 = sim.output.figure_axe() -title = ('intermittency, solver '+sim.output.name_solver+ -f', nh = {sim.param.nx:5d}'+ -', c = {:.4g}, f = {:.4g}'.format(np.sqrt(sim.param.c2), sim.param.f) +title = ( + "intermittency, solver " + + sim.output.name_solver + + f", nh = {sim.param.nx:5d}" + + ", c = {:.4g}, f = {:.4g}".format(np.sqrt(sim.param.c2), sim.param.f) ) ax1.set_title(title) ax1.hold(True) -ax1.set_xscale('linear') -ax1.set_yscale('linear') +ax1.set_xscale("linear") +ax1.set_yscale("linear") -ax1.set_xlabel('order q') -ax1.set_ylabel(r'$\zeta_q$') +ax1.set_xlabel("order q") +ax1.set_ylabel(r"$\zeta_q$") ax1.plot(orders, expos) -ax1.plot(orders, expos_K41, 'k--') +ax1.plot(orders, expos_K41, "k--") solveq2d.show() diff --git a/scripts/plot_results/plot_many_runs.py b/scripts/plot_results/plot_many_runs.py index 2193264c..7fd47c67 100644 --- a/scripts/plot_results/plot_many_runs.py +++ b/scripts/plot_results/plot_many_runs.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -#coding=utf8 +# coding=utf8 import numpy as np import glob @@ -14,63 +14,55 @@ # dir_base = 'Pure_standing_waves_512x512' dir_base = ( -# '/scratch/augier/' -'/home/pierre/' -'Results_for_article_SW1l/' -'Pure_standing_waves_1024x1024' + # '/scratch/augier/' + "/home/pierre/" + "Results_for_article_SW1l/" + "Pure_standing_waves_1024x1024" ) - - - - - - set_of_dir_results = solveq2d.SetOfDirResults(dir_base=dir_base) values_c = np.array([10, 1000]) # values_c2 = set_of_dir_results.values_c2 -values_c2 = values_c**2 +values_c2 = values_c ** 2 values_solver = set_of_dir_results.values_solver # values_solver = ['SW1l'] tmin = 80 -tuple_loop = [(c2,name_solver) - for c2 in values_c2 - for name_solver in values_solver] +tuple_loop = [ + (c2, name_solver) for c2 in values_c2 for name_solver in values_solver +] for c2, name_solver in tuple_loop: path_dir = set_of_dir_results.one_path_from_values( - solver=name_solver, - c2=c2, - FORCING=True) + solver=name_solver, c2=c2, FORCING=True + ) if path_dir is None: continue sim = solveq2d.load_state_phys_file(t_approx=2620, name_dir=path_dir) - sim.output.spect_energy_budg.plot(tmin=tmin, tmax=500, delta_t=0.) - sim.output.spectra.plot2D(tmin=tmin, tmax=500, delta_t=0., - coef_compensate=5./3) + sim.output.spect_energy_budg.plot(tmin=tmin, tmax=500, delta_t=0.0) + sim.output.spectra.plot2D( + tmin=tmin, tmax=500, delta_t=0.0, coef_compensate=5.0 / 3 + ) # sim.output.spatial_means.plot() sim.output.prob_dens_func.plot(tmin=tmin) - sim.output.increments.plot(tmin=tmin, tmax=None, delta_t=0., - order=6, yscale='log') + sim.output.increments.plot( + tmin=tmin, tmax=None, delta_t=0.0, order=6, yscale="log" + ) - # sim.output.increments.plot_pdf(tmin=tmin, tmax=160.25,key_var='ux', + # sim.output.increments.plot_pdf(tmin=tmin, tmax=160.25,key_var='ux', # order=4) sim.output.increments.plot_Kolmo(tmin=tmin) - - - # sim.output.phys_fields.plot(key_field='rot') # sim.output.phys_fields.plot(key_field='eta') diff --git a/scripts/plot_results/plot_many_things.py b/scripts/plot_results/plot_many_things.py index 5237ab59..ebf0035f 100644 --- a/scripts/plot_results/plot_many_things.py +++ b/scripts/plot_results/plot_many_things.py @@ -1,16 +1,13 @@ - - import matplotlib.pylab as plt from solveq2d import solveq2d - name_dir = ( -'/scratch/augier/' -'/Results_SW1lw' -'/Pure_standing_waves_7680x7680' -'/SE2D_SW1lwaves_forcingw_L=50.x50._7680x7680_c=40_f=0_2013-08-06_12-26-05' + "/scratch/augier/" + "/Results_SW1lw" + "/Pure_standing_waves_7680x7680" + "/SE2D_SW1lwaves_forcingw_L=50.x50._7680x7680_c=40_f=0_2013-08-06_12-26-05" ) tmin = 194 @@ -22,21 +19,20 @@ # sim.output.phys_fields.plot(key_field='eta') - sim = solveq2d.create_sim_plot_from_dir(name_dir=name_dir) # sim.output.print_stdout.plot() # sim.output.spatial_means.plot() -# sim.output.spectra.plot1D(tmin=tmin, tmax=tmax, delta_t=0., +# sim.output.spectra.plot1D(tmin=tmin, tmax=tmax, delta_t=0., # coef_compensate=2.) -# sim.output.spectra.plot2D(tmin=tmin, tmax=tmax, delta_t=0., +# sim.output.spectra.plot2D(tmin=tmin, tmax=tmax, delta_t=0., # coef_compensate=2.) # sim.output.spect_energy_budg.plot(tmin=tmin, tmax=tmax, delta_t=0.) -# sim.output.increments.plot(tmin=tmin, tmax=None, delta_t=0., +# sim.output.increments.plot(tmin=tmin, tmax=None, delta_t=0., # order=4, yscale='log') # sim.output.prob_dens_func.plot(tmin=tmin) @@ -45,5 +41,3 @@ solveq2d.show() - - diff --git a/scripts/plot_results/plot_profil_phys.py b/scripts/plot_results/plot_profil_phys.py index 5f89eb79..4c05f896 100644 --- a/scripts/plot_results/plot_profil_phys.py +++ b/scripts/plot_results/plot_profil_phys.py @@ -1,40 +1,41 @@ - import numpy as np import matplotlib.pylab as plt from solveq2d import solveq2d + def jump_var_delta(sim, var, deltar): nx = sim.param.nx - dx = sim.param.Lx/nx - idelta = int(round(deltar/dx)) + dx = sim.param.Lx / nx + idelta = int(round(deltar / dx)) jump = -var # jump = np.zeros(var.shape) - jump[:,0:nx-idelta] += var[:,idelta:nx] - jump[:,nx-idelta:nx] += var[:,0:idelta] + jump[:, 0 : nx - idelta] += var[:, idelta:nx] + jump[:, nx - idelta : nx] += var[:, 0:idelta] return jump + name_dir = ( -'/scratch/augier' -# '/home/pierre' -'/Results_SW1lw' -'/Pure_standing_waves_512x512/' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._512x512_c=10_f=0_2013-06-12_00-02-04' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._512x512_c=200_f=0_2013-06-12_00-39-17' -'SE2D_SW1lwaves_forcingw_L=50.x50._512x512_c=700_f=0_2013-06-12_12-40-31' -# '/Pure_standing_waves_1024x1024/' -# '/Pure_standing_waves_512x512/' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=10_f=0_2013-06-12_19-34-30' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=20_f=0_2013-06-12_20-33-19' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=50_f=0_2013-06-12_21-35-00' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=100_f=0_2013-06-12_22-44-21' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=200_f=0_2013-06-13_00-33-41' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=400_f=0_2013-06-13_03-42-31' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=700_f=0_2013-06-13_09-32-03' -# 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=1000_f=0_2013-06-13_19-31-10' + "/scratch/augier" + # '/home/pierre' + "/Results_SW1lw" + "/Pure_standing_waves_512x512/" + # 'SE2D_SW1lwaves_forcingw_L=50.x50._512x512_c=10_f=0_2013-06-12_00-02-04' + # 'SE2D_SW1lwaves_forcingw_L=50.x50._512x512_c=200_f=0_2013-06-12_00-39-17' + "SE2D_SW1lwaves_forcingw_L=50.x50._512x512_c=700_f=0_2013-06-12_12-40-31" + # '/Pure_standing_waves_1024x1024/' + # '/Pure_standing_waves_512x512/' + # 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=10_f=0_2013-06-12_19-34-30' + # 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=20_f=0_2013-06-12_20-33-19' + # 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=50_f=0_2013-06-12_21-35-00' + # 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=100_f=0_2013-06-12_22-44-21' + # 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=200_f=0_2013-06-13_00-33-41' + # 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=400_f=0_2013-06-13_03-42-31' + # 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=700_f=0_2013-06-13_09-32-03' + # 'SE2D_SW1lwaves_forcingw_L=50.x50._1024x1024_c=1000_f=0_2013-06-13_19-31-10' ) -deltar = 1. +deltar = 1.0 divJ_limit = 1 sim = solveq2d.load_state_phys_file(t_approx=262, name_dir=name_dir) @@ -46,20 +47,19 @@ def jump_var_delta(sim, var, deltar): name_solver = sim.output.name_solver - # sim.output.phys_fields.plot(key_field='eta') # sim.output.phys_fields.plot(key_field='div') # sim.output.phys_fields.plot(key_field='Jx') # sim.output.phys_fields.plot(key_field='Jy') -eta = sim.state.state_phys['eta'] -h = eta+1 +eta = sim.state.state_phys["eta"] +h = eta + 1 -ux = sim.state.state_phys['ux'] -uy = sim.state.state_phys['uy'] +ux = sim.state.state_phys["ux"] +uy = sim.state.state_phys["uy"] -Jx = h*ux -Jy = h*uy +Jx = h * ux +Jy = h * uy oper = sim.oper @@ -76,50 +76,45 @@ def jump_var_delta(sim, var, deltar): px_h = oper.ifft2(px_h_fft) py_h = oper.ifft2(py_h_fft) -gradh2 = px_h**2 + py_h**2 +gradh2 = px_h ** 2 + py_h ** 2 # sim.output.phys_fields.plot(field=gradh2) -gradh2_limit = 60./c2 +gradh2_limit = 60.0 / c2 -Usx = divJ*px_h/gradh2 -Usy = divJ*py_h/gradh2 +Usx = divJ * px_h / gradh2 +Usy = divJ * py_h / gradh2 -Usx[abs(gradh2)maxv] = maxv - var[var maxv] = maxv + var[var < minv] = minv + maxMach = 1.5 minMach = 0.3 -limit_maxmin(Us_over_c, 0., 2.) - -sim.output.phys_fields.plot(field=Us_over_c, key_field='Us/c', - nb_contours=20) +limit_maxmin(Us_over_c, 0.0, 2.0) +sim.output.phys_fields.plot(field=Us_over_c, key_field="Us/c", nb_contours=20) - - -Machx2 = (ux-Usx)**2/c2 -Machy2 = (uy-Usy)**2/c2 +Machx2 = (ux - Usx) ** 2 / c2 +Machy2 = (uy - Usy) ** 2 / c2 Mach = np.sqrt(Machx2 + Machy2) -Mach[abs(gradh2) phi0), u1, u0) ug = ux_from_lats(lats1r) ug.shape = (oper.nlat, 1) -ux = ug * oper.create_array_spat(value=1.) # broadcast to shape (nlats,nlons) +ux = ug * oper.create_array_spat(value=1.0) # broadcast to shape (nlats,nlons) # uy = self.state_phys.get_var("uy") def integrand_gh(phi, a, omega): @@ -124,17 +124,17 @@ def h_from_eta(eta): etag = eta_from_h(h) etag.shape = (oper.nlat, 1) -eta = etag * oper.create_array_spat(value=1.) # broadcast to shape (nlats,nlons) +eta = etag * oper.create_array_spat(value=1.0) # broadcast to shape (nlats,nlons) # Height perturbation. -alpha = 1. / 3. -beta = 1. / 15. -hamp = 120. # amplitude of height perturbation to zonal jet +alpha = 1.0 / 3.0 +beta = 1.0 / 15.0 +hamp = 120.0 # amplitude of height perturbation to zonal jet hbump = ( hamp * np.cos(LATS) - * np.exp(-(LONS / alpha) ** 2) - * np.exp(-((phi2 - LATS) / beta) ** 2) + * np.exp(-((LONS / alpha) ** 2)) + * np.exp(-(((phi2 - LATS) / beta) ** 2)) ) eta += g * hbump / params.c2 diff --git a/scripts/util/modif_resol_all_dir.py b/scripts/util/modif_resol_all_dir.py index 907faa9a..06cdc234 100644 --- a/scripts/util/modif_resol_all_dir.py +++ b/scripts/util/modif_resol_all_dir.py @@ -3,11 +3,11 @@ from solveq2d import solveq2d t_approx = None -coef_modif_resol=2 -dir_base = '~/Storage/Results_SW1lw/Pure_standing_waves_1920x1920' +coef_modif_resol = 2 +dir_base = "~/Storage/Results_SW1lw/Pure_standing_waves_1920x1920" -print('run the function solveq2d.modif_resolution_all_dir') +print("run the function solveq2d.modif_resolution_all_dir") -solveq2d.modif_resolution_all_dir(t_approx=t_approx, - coef_modif_resol=coef_modif_resol, - dir_base=dir_base) +solveq2d.modif_resolution_all_dir( + t_approx=t_approx, coef_modif_resol=coef_modif_resol, dir_base=dir_base +) diff --git a/scripts/util/modif_resolution.py b/scripts/util/modif_resolution.py index 30ccc8dd..566ca153 100644 --- a/scripts/util/modif_resolution.py +++ b/scripts/util/modif_resolution.py @@ -1,21 +1,18 @@ #!/usr/bin/env python -#coding=utf8 -# +# coding=utf8 +# # run modif_resolution.py # python modif_resolution.py from solveq2d import solveq2d name_dir = ( - '~/Storage/Results_SW1lw/Pure_standing_waves_3840x3840/'+ - 'SE2D_SW1lwaves_forcingw_L=50.x50._3840x3840_c=40_f=0_2013-07-23_15-47-25' + "~/Storage/Results_SW1lw/Pure_standing_waves_3840x3840/" + + "SE2D_SW1lwaves_forcingw_L=50.x50._3840x3840_c=40_f=0_2013-07-23_15-47-25" ) -t_approx = 10.e10 -coef_modif_resol=2 +t_approx = 10.0e10 +coef_modif_resol = 2 solveq2d.modif_resolution_from_dir( - name_dir, t_approx=t_approx, - coef_modif_resol=coef_modif_resol, - PLOT=False) - - + name_dir, t_approx=t_approx, coef_modif_resol=coef_modif_resol, PLOT=False +) diff --git a/scripts/util/resume_from_path.py b/scripts/util/resume_from_path.py index 9668f761..09140385 100644 --- a/scripts/util/resume_from_path.py +++ b/scripts/util/resume_from_path.py @@ -1,9 +1,10 @@ import argparse import fluidsim as fls -parser = argparse.ArgumentParser(description='Resume a Fluidsim simulation from a path') -parser.add_argument('path', - help='path to an incomplete simulation directory') +parser = argparse.ArgumentParser( + description="Resume a Fluidsim simulation from a path" +) +parser.add_argument("path", help="path to an incomplete simulation directory") args = parser.parse_args() sim = fls.load_state_phys_file(args.path, modif_save_params=False)