Clean up examples

hg hash: a3fc57f974d248125c0dd7b478f99bd3be536d30

Author: kburns

examples/bvp/1d_lane_emden/lane_emden.py

@@ -52,7 +52,6 @@
 
 # Setup problem
 problem = de.NLBVP(domain, variables=['f', 'fx', 'R'], ncc_cutoff=ncc_cutoff)
-problem.meta[:]['x']['dirichlet'] = True
 problem.parameters['n'] = n
 problem.add_equation("x*dx(fx) + 2*fx = -x*(R**2)*(f**n)", tau=False)
 problem.add_equation("fx - dx(f) = 0")

examples/evp/1d_rayleigh_benard/rayleigh_benard.py

@@ -34,7 +34,6 @@
 # 2D Boussinesq hydrodynamics, with no-slip boundary conditions
 # Use substitutions for x and t derivatives
 problem = de.EVP(domain, variables=['p','b','u','w','bz','uz','wz'], eigenvalue='omega')
-problem.meta[:]['z']['dirichlet'] = True
 problem.parameters['P'] = (Rayleigh * Prandtl)**(-1/2)
 problem.parameters['R'] = (Rayleigh / Prandtl)**(-1/2)
 problem.parameters['F'] = F = 1

examples/evp/1d_waves_on_string/waves_on_a_string.py

@@ -29,7 +29,6 @@
 
 # Problem
 problem = de.EVP(domain, variables=['u', 'u_x'],eigenvalue='l')
-problem.meta[:]['x']['dirichlet'] = True
 problem.add_equation("l*u + dx(u_x) = 0")
 problem.add_equation("u_x - dx(u) = 0")
 problem.add_bc("left(u) = 0")

examples/ivp/1d_kdv_burgers/kdv_burgers.py

@@ -45,7 +45,7 @@
 ux.differentiate(0, out=uxx)
 
 # Store data for final plot
-u.set_scales(1, keep_data=True)
+u.set_scales(1)
 u_list = [np.copy(u['g'])]
 t_list = [solver.sim_time]
 
@@ -54,7 +54,7 @@
 while solver.ok:
     solver.step(dt)
     if solver.iteration % 20 == 0:
-        u.set_scales(1, keep_data=True)
+        u.set_scales(1)
         u_list.append(np.copy(u['g']))
         t_list.append(solver.sim_time)
     if solver.iteration % 100 == 0:

examples/ivp/2d_boundary_forced_waves/boundary_forced_waves.py

@@ -1,5 +1,4 @@
-"""boundary_forced_waves.py
-
+"""
 Solve the 2D linear wave equation forced on the right y-boundary by an
 x and t-dependent forcing function:
 
@@ -19,6 +18,7 @@
 with n = 10, Lf = 1, and freq = 10.
 
 """
+
 import time
 import numpy as np
 import dedalus.public as de
@@ -27,81 +27,78 @@
 import logging
 logger = logging.getLogger(__name__)
 
+
+# Parameters
 nx = 100
 ny = 100
-c = 1.
+c = 1.0
+dt = 0.01
 
-x = de.Fourier('x',nx,interval=[0.,2*np.pi])
-y = de.Chebyshev('y',ny,interval=[0.,2*np.pi])
-domain = de.Domain([x,y],grid_dtype='float')
+# Bases and domain
+x_basis = de.Fourier('x', nx, interval=(0, 2*np.pi))
+y_basis = de.Chebyshev('y', ny, interval=(0, 2*np.pi))
+domain = de.Domain([x_basis, y_basis], grid_dtype='float')
 
 def BoundaryForcing(*args):
-    """this function applies its arguments and returns the forcing
-
-    """
+    """This function applies its arguments and returns the forcing"""
     t = args[0].value # this is a scalar; we use .value to get its value
     x = args[1].data # this is an array; we use .data to get its values
     ampl = args[2].value
     freq = args[3].value
-
     return ampl*cos_bump(x)*np.cos(t*freq)
 
 def cos_bump(x, Lx=2*np.pi, n=10):
-    """a simple, smooth bump function.
-
-    """
+    """A simple, smooth bump function."""
     return (((1-np.cos(2*np.pi/Lx * x))/2)**n)
 
-
 def Forcing(*args, domain=domain, F=BoundaryForcing):
-    """This function takes arguments *args, a function F, and a domain and
+    """
+    This function takes arguments *args, a function F, and a domain and
     returns a Dedalus GeneralFunction that can be applied.
-
     """
     return de.operators.GeneralFunction(domain, layout='g', func=F, args=args)
 
-# now we make it parseable, so the symbol BF can be used in equations
+# Now we make it parseable, so the symbol BF can be used in equations
 # and the parser will know what it is.
 de.operators.parseables['BF'] = Forcing
 
-waves = de.IVP(domain,variables=['f','f_t','f_y'])
+# Problem
+waves = de.IVP(domain, variables=['f','f_t','f_y'])
 waves.parameters['csq'] = c**2
 waves.parameters['ampl'] = 1.
 waves.parameters['freq'] = 10.
 waves.add_equation("dt(f_t) - csq*(dx(dx(f)) + dy(f_y)) = 0")
 waves.add_equation("f_t - dt(f) = 0")
 waves.add_equation("f_y - dy(f) = 0")
-
 # note that we apply the right() operator to the forcing
 waves.add_bc("left(f) = 0")
 waves.add_bc("right(f) = right(BF(t,x,ampl,freq))")
 
+# Solver
 solver = waves.build_solver(de.timesteppers.RK443)
-
 solver.stop_sim_time = 10.
-solver.stop_wall_time = 1000.
-solver.stop_iteration = np.inf
-dt = 0.01
 
-# analysis
+# Analysis
 analysis_tasks = []
 check = solver.evaluator.add_file_handler('checkpoints', iter=10, max_writes=200)
 check.add_system(solver.state)
 analysis_tasks.append(check)
 
+# Main loop
 start_time = time.time()
 while solver.ok:
     solver.step(dt)
     if (solver.iteration-1) % 10 == 0:
         logger.info("Time step {}".format(solver.iteration))
 end_time = time.time()
+
 # Print statistics
 logger.info('Total time: %f' %(end_time-start_time))
 logger.info('Iterations: %i' %solver.iteration)
 logger.info('Average timestep: %f' %(solver.sim_time/solver.iteration))
 
+# Merge output
 logger.info('beginning join operation')
-
 for task in analysis_tasks:
     logger.info(task.base_path)
     post.merge_analysis(task.base_path)

examples/ivp/2d_rayleigh_benard/merge.py

@@ -2,7 +2,7 @@
 Plot planes from joint analysis files.
 
 Usage:
-    plot_2d_series.py <files>... [--output=<dir>]
+    plot_slices.py <files>... [--output=<dir>]
 
 Options:
     --output=<dir>  Output directory [default: ./frames]

…ivp/2d_rayleigh_benard/plot_2d_series.py …es/ivp/2d_rayleigh_benard/plot_slices.pyexamples/ivp/2d_rayleigh_benard/plot_2d_series.py renamed to examples/ivp/2d_rayleigh_benard/plot_slices.py examples/ivp/2d_rayleigh_benard/plot_2d_series.py renamed to examples/ivp/2d_rayleigh_benard/plot_slices.py

@@ -5,22 +5,34 @@
 conditions.  The equations are scaled in units of the buoyancy time (Fr = 1).
 
 This script can be ran serially or in parallel, and uses the built-in analysis
-framework to save data snapshots in HDF5 files.  The `merge.py` script in this
-folder can be used to merge distributed analysis sets from parallel runs,
-and the `plot_2d_series.py` script can be used to plot the snapshots.
+framework to save data snapshots in HDF5 files.  The `merge_procs` command can
+be used to merge distributed analysis sets from parallel runs, and the
+`plot_slices.py` script can be used to plot the snapshots.
 
 To run, merge, and plot using 4 processes, for instance, you could use:
     $ mpiexec -n 4 python3 rayleigh_benard.py
-    $ mpiexec -n 4 python3 merge.py snapshots
-    $ mpiexec -n 4 python3 plot_2d_series.py snapshots/*.h5
+    $ mpiexec -n 4 python3 -m dedalus merge_procs snapshots
+    $ mpiexec -n 4 python3 plot_slices.py snapshots/*.h5
 
-The simulation should take a few process-minutes to run.
+This script can restart the simulation from the last save of the original
+output to extend the integration.  This requires that the output files from
+the original simulation are merged, and the last is symlinked or copied to
+`restart.h5`.
+
+To run the original example and the restart, you could use:
+    $ mpiexec -n 4 python3 rayleigh_benard.py
+    $ mpiexec -n 4 python3 -m dedalus merge_procs snapshots
+    $ ln -s snapshots/snapshots_s2.h5 restart.h5
+    $ mpiexec -n 4 python3 rayleigh_benard.py
+
+The simulations should take a few process-minutes to run.
 
 """
 
 import numpy as np
 from mpi4py import MPI
 import time
+import pathlib
 
 from dedalus import public as de
 from dedalus.extras import flow_tools
@@ -64,34 +76,46 @@
 solver = problem.build_solver(de.timesteppers.RK222)
 logger.info('Solver built')
 
-# Initial conditions
-x = domain.grid(0)
-z = domain.grid(1)
-b = solver.state['b']
-bz = solver.state['bz']
-
-# Random perturbations, initialized globally for same results in parallel
-gshape = domain.dist.grid_layout.global_shape(scales=1)
-slices = domain.dist.grid_layout.slices(scales=1)
-rand = np.random.RandomState(seed=42)
-noise = rand.standard_normal(gshape)[slices]
-
-# Linear background + perturbations damped at walls
-zb, zt = z_basis.interval
-pert =  1e-3 * noise * (zt - z) * (z - zb)
-b['g'] = F * pert
-b.differentiate('z', out=bz)
-
-# Initial timestep
-dt = 0.125
+# Initial conditions or restart
+if not pathlib.Path('restart.h5').exists():
+
+    # Initial conditions
+    x = domain.grid(0)
+    z = domain.grid(1)
+    b = solver.state['b']
+    bz = solver.state['bz']
+
+    # Random perturbations, initialized globally for same results in parallel
+    gshape = domain.dist.grid_layout.global_shape(scales=1)
+    slices = domain.dist.grid_layout.slices(scales=1)
+    rand = np.random.RandomState(seed=42)
+    noise = rand.standard_normal(gshape)[slices]
+
+    # Linear background + perturbations damped at walls
+    zb, zt = z_basis.interval
+    pert =  1e-3 * noise * (zt - z) * (z - zb)
+    b['g'] = F * pert
+    b.differentiate('z', out=bz)
+
+    # Timestepping and output
+    dt = 0.125
+    stop_sim_time = 25
+    fh_mode = 'overwrite'
+
+else:
+    # Restart
+    write, last_dt = solver.load_state('restart.h5', -1)
+
+    # Timestepping and output
+    dt = last_dt
+    stop_sim_time = 50
+    fh_mode = 'append'
 
 # Integration parameters
-solver.stop_sim_time = 25
-solver.stop_wall_time = 30 * 60.
-solver.stop_iteration = np.inf
+solver.stop_sim_time = stop_sim_time
 
 # Analysis
-snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.25, max_writes=50)
+snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.25, max_writes=50, mode=fh_mode)
 snapshots.add_system(solver.state)
 
 # CFL