Read about our upcoming Code of Conduct on this issue

Commit 015b5cf5 authored by Ashwin V. Mohanan's avatar Ashwin V. Mohanan
Browse files

Improve diagnostics, set oddball mode to zero in Burgers

parent 166cbf9bbdfd
Pipeline #22795 passed with stage
in 16 minutes and 47 seconds
...@@ -8,6 +8,13 @@ Notes ...@@ -8,6 +8,13 @@ Notes
- This contrast is less evident if nx=64, 128 etc. which is indeed puzzling. - This contrast is less evident if nx=64, 128 etc. which is indeed puzzling.
- As suggested in the reference below the foolproof solution is to apply - As suggested in the reference below the foolproof solution is to apply
dealiasing which guarantees energy conservation. dealiasing which guarantees energy conservation.
- The time scheme also plays a role in influencing the final energy, but is
less (?) likely to affect aliasing errors directly.
- It seems that the oddball mode was not set to zero in FluidSim, but this is
perhaps not a big problem since we always use some dealiasing?
Examples
--------
.. seealso: https://kth-nek5000.github.io/kthNekBook/_notebooks/burgers.html .. seealso: https://kth-nek5000.github.io/kthNekBook/_notebooks/burgers.html
...@@ -17,7 +24,7 @@ import numpy as np ...@@ -17,7 +24,7 @@ import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
def solve(Simul, nx, coef_dealiasing): def solve(Simul, nx, time_scheme, coef_dealiasing):
params = Simul.create_default_params() params = Simul.create_default_params()
params.output.sub_directory = "bench_skew_sym" params.output.sub_directory = "bench_skew_sym"
...@@ -29,8 +36,8 @@ def solve(Simul, nx, coef_dealiasing): ...@@ -29,8 +36,8 @@ def solve(Simul, nx, coef_dealiasing):
params.nu_2 = 0.0 params.nu_2 = 0.0
params.time_stepping.USE_CFL = False params.time_stepping.USE_CFL = False
params.time_stepping.deltat0 = 1e-2 params.time_stepping.deltat0 = 1e-2
params.time_stepping.t_end = 1.0 params.time_stepping.t_end = 1.1
params.time_stepping.type_time_scheme = "RK2" params.time_stepping.type_time_scheme = time_scheme
params.init_fields.type = "in_script" params.init_fields.type = "in_script"
params.output.periods_print.print_stdout = 0.5 params.output.periods_print.print_stdout = 0.5
...@@ -43,33 +50,62 @@ def solve(Simul, nx, coef_dealiasing): ...@@ -43,33 +50,62 @@ def solve(Simul, nx, coef_dealiasing):
u = np.sin(x) u = np.sin(x)
u_fft = sim.oper.fft(u) u_fft = sim.oper.fft(u)
# Set "oddball mode" to zero
u_fft[sim.oper.nkx-1] = 0.
sim.state.init_statephys_from(u=u) sim.state.init_statephys_from(u=u)
sim.state.init_statespect_from(u_fft=u_fft) sim.state.init_statespect_from(u_fft=u_fft)
E_initial = sim.output.compute_energy()
# Plot once # Plot once
sim.output.init_with_initialized_state() sim.output.init_with_initialized_state()
sim.output.phys_fields.plot(field="u") sim.output.phys_fields.plot(field="u")
ax = plt.gca() ax = plt.gca()
ax.grid()
ax.set(xlim=(0, Lx))
sim.time_stepping.start() sim.time_stepping.start()
E_final = sim.output.compute_energy()
print("Final Energy / Initial Energy =", E_final / E_initial)
ax.plot(x, sim.state.state_phys.get_var("u")) ax.plot(x, sim.state.state_phys.get_var("u"))
plt.show() plt.show()
@click.command() @click.command()
@click.option("--solver", type=click.Choice(["conv", "skew_sym"])) @click.option("--solver", type=click.Choice(["conv", "skew_sym"]))
@click.option("--nx", type=int, default=101) @click.option("--nx", type=int, default=128)
@click.option("--ts", type=click.Choice(["Euler", "RK2", "RK4"]), default="Euler")
@click.option("--dealias/--no-dealias", type=bool, default=False) @click.option("--dealias/--no-dealias", type=bool, default=False)
def run(solver, nx, dealias): def run(solver, nx, ts, dealias):
"""Execute Burgers solvers. Compare the following runs:
\b
# No dealiasing
>>> ./run.py --nx 128 --solver conv
>>> ./run.py --nx 128 --solver skew_sym
\b
# With dealiasing
>>> ./run.py --nx 192 --dealias --solver conv
>>> ./run.py --nx 192 --dealias --solver skew_sym
"""
if solver == "conv": if solver == "conv":
from fluidsim.solvers.burgers1d.solver import Simul from fluidsim.solvers.burgers1d.solver import Simul
else: else:
from fluidsim.solvers.burgers1d.skew_sym.solver import Simul from fluidsim.solvers.burgers1d.skew_sym.solver import Simul
solve(Simul, nx, 2./3 if dealias else 1.0) solve(Simul, nx, ts, 2. / 3 if dealias else 1.0)
if __name__ == "__main__": if __name__ == "__main__":
params = {
"figure.figsize": (14, 5),
"grid.alpha": 0.6,
}
plt.rcParams.update(params)
run() run()
...@@ -58,4 +58,6 @@ class Simul(SimulBurgers): ...@@ -58,4 +58,6 @@ class Simul(SimulBurgers):
f_fft[:] += f_fft_cons f_fft[:] += f_fft_cons
self.oper.dealiasing(f_fft) self.oper.dealiasing(f_fft)
# Set "oddball mode" to zero
f_fft[self.oper.nkx-1] = 0.
return tendencies_fft return tendencies_fft
...@@ -89,6 +89,8 @@ class Simul(SimulBasePseudoSpectral): ...@@ -89,6 +89,8 @@ class Simul(SimulBasePseudoSpectral):
f_fft = tendencies_fft.get_var("u_fft") f_fft = tendencies_fft.get_var("u_fft")
self.oper.fft_as_arg(f_signal, f_fft) self.oper.fft_as_arg(f_signal, f_fft)
self.oper.dealiasing(f_fft) self.oper.dealiasing(f_fft)
# Set "oddball mode" to zero
f_fft[self.oper.nkx-1] = 0.
return tendencies_fft return tendencies_fft
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment