-
Notifications
You must be signed in to change notification settings - Fork 33
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
2 changed files
with
121 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,118 @@ | ||
#!/usr/bin/env python3 | ||
# -*- coding: utf-8 -*- | ||
""" | ||
Dedalus script simulating 2D horizontally-periodic Rayleigh-Benard convection | ||
using Spectral Deferred Corrections. | ||
This script demonstrates solving a 2D Cartesian initial value problem. It can | ||
be ran serially or in parallel, and uses the built-in analysis framework to save | ||
data snapshots to HDF5 files. The `plot_snapshots.py` script can be used to | ||
produce plots from the saved data. It should take about 5 cpu-minutes to run. | ||
The problem is non-dimensionalized using the box height and freefall time, so | ||
the resulting thermal diffusivity and viscosity are related to the Prandtl | ||
and Rayleigh numbers as: | ||
kappa = (Rayleigh * Prandtl)**(-1/2) | ||
nu = (Rayleigh / Prandtl)**(-1/2) | ||
For incompressible hydro with two boundaries, we need two tau terms for each the | ||
velocity and buoyancy. Here we choose to use a first-order formulation, putting | ||
one tau term each on auxiliary first-order gradient variables and the others in | ||
the PDE, and lifting them all to the first derivative basis. This formulation puts | ||
a tau term in the divergence constraint, as required for this geometry. | ||
To run using e.g. 4 processes: | ||
$ mpiexec -n 4 python rayleighBenardSDC.py | ||
Then use https://github.com/DedalusProject/dedalus/blob/master/examples/ivp_2d_rayleigh_benard/plot_snapshots.py | ||
for plotting. | ||
""" | ||
import numpy as np | ||
import dedalus.public as d3 | ||
from pySDC.playgrounds.dedalus.sdc import SpectralDeferredCorrectionIMEX | ||
import logging | ||
logger = logging.getLogger(__name__) | ||
|
||
useSDC = True | ||
SpectralDeferredCorrectionIMEX.setParameters( | ||
nSweeps=4, | ||
nNodes=4, | ||
implSweep="MIN-SR-FLEX", | ||
explSweep="PIC") | ||
|
||
# Parameters | ||
Lx, Lz = 4, 1 | ||
Nx, Nz = 512, 128 | ||
Rayleigh = 2e6 | ||
Prandtl = 1 | ||
dealias = 3/2 | ||
stop_sim_time = 50 | ||
timestepper = SpectralDeferredCorrectionIMEX if useSDC else d3.RK222 | ||
timestep = 1e-2/4 | ||
dtype = np.float64 | ||
|
||
# Bases | ||
coords = d3.CartesianCoordinates('x', 'z') | ||
dist = d3.Distributor(coords, dtype=dtype) | ||
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx), dealias=dealias) | ||
zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz), dealias=dealias) | ||
|
||
# Fields | ||
p = dist.Field(name='p', bases=(xbasis,zbasis)) | ||
b = dist.Field(name='b', bases=(xbasis,zbasis)) | ||
u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis)) | ||
tau_p = dist.Field(name='tau_p') | ||
tau_b1 = dist.Field(name='tau_b1', bases=xbasis) | ||
tau_b2 = dist.Field(name='tau_b2', bases=xbasis) | ||
tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis) | ||
tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis) | ||
|
||
# Substitutions | ||
kappa = (Rayleigh * Prandtl)**(-1/2) | ||
nu = (Rayleigh / Prandtl)**(-1/2) | ||
x, z = dist.local_grids(xbasis, zbasis) | ||
ex, ez = coords.unit_vector_fields(dist) | ||
lift_basis = zbasis.derivative_basis(1) | ||
lift = lambda A: d3.Lift(A, lift_basis, -1) | ||
grad_u = d3.grad(u) + ez*lift(tau_u1) # First-order reduction | ||
grad_b = d3.grad(b) + ez*lift(tau_b1) # First-order reduction | ||
|
||
# Problem | ||
# First-order form: "div(f)" becomes "trace(grad_f)" | ||
# First-order form: "lap(f)" becomes "div(grad_f)" | ||
problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals()) | ||
problem.add_equation("trace(grad_u) + tau_p = 0") | ||
problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - u@grad(b)") | ||
problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = - u@grad(u)") | ||
problem.add_equation("b(z=0) = Lz") | ||
problem.add_equation("u(z=0) = 0") | ||
problem.add_equation("b(z=Lz) = 0") | ||
problem.add_equation("u(z=Lz) = 0") | ||
problem.add_equation("integ(p) = 0") # Pressure gauge | ||
|
||
# Solver | ||
solver = problem.build_solver(timestepper) | ||
solver.stop_sim_time = stop_sim_time | ||
|
||
# Initial conditions | ||
b.fill_random('g', seed=42, distribution='normal', scale=1e-3) # Random noise | ||
b['g'] *= z * (Lz - z) # Damp noise at walls | ||
b['g'] += Lz - z # Add linear background | ||
|
||
# Analysis | ||
snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.25, max_writes=50) | ||
snapshots.add_task(b, name='buoyancy') | ||
snapshots.add_task(-d3.div(d3.skew(u)), name='vorticity') | ||
|
||
# Main loop | ||
try: | ||
logger.info('Starting main loop') | ||
while solver.proceed: | ||
solver.step(timestep) | ||
if (solver.iteration-1) % 10 == 0: | ||
logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep)) | ||
except: | ||
logger.error('Exception raised, triggering end of main loop.') | ||
raise | ||
finally: | ||
solver.log_stats() |