How to implement lateral thickness variations in 2D Rayleigh–Benard convection example
saurabh shukla
Hi all,
I am trying to modify the ivp_2d_rayleigh_benard.py example in Dedalus to include lateral variations in the thickness of the layer.
The approach my advisor suggested is to map the vertical coordinate to a constant thickness one:
where is the depth and is the thickness of the ice/fluid layer as a function of the horizontal coordinate .
The idea is that these rescaled equations could then be implemented in Dedalus, starting from the Rayleigh–Benard convection setup.
My questions are:
What is the best way to implement this coordinate transformation in Dedalus v3?
Should I explicitly rewrite the PDEs using the chain rule with , or
Can Dedalus's coordinate system / metric framework handle this transformation more naturally?
Has anyone already implemented lateral thickness variations (or a similar coordinate scaling) in Dedalus that I could look at as a reference?
At the moment, I'm stuck on how to correctly define the derivatives and operators after this coordinate transformation, so any guidance or example code would be really helpful.
Thanks!!!
Saurabh
PS : Here is the snippet of my code (modified version of the 2D rayleigh-benard code)
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>..
"""
Dedalus script simulating 2D Rayleigh-Benard convection with variable thickness.
Modified from the original 2D RB example to implement coordinate transformation
Z = z/d(x) as suggested by professor, where d(x) is the variable thickness function.
This maps the variable-thickness physical domain to a constant computational domain.
Key changes from original (I made these changes to fit in the model):
1. Changed vertical coordinate from z to Z (computational coordinate)
2. Added thickness function d(x)
3. Implemented coordinate transformation in all derivatives
4. Modified boundary conditions for variable thickness
"""
import numpy as np
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)
# Parameters
Lx = 4
Nx, Nz = 256, 64
Rayleigh = 2e6
Prandtl = 1
dealias = 3/2
stop_sim_time = 50
timestepper = d3.RK222
max_timestep = 0.125
dtype = np.float64
# ADDED: Variable thickness function
def thickness(x):
"""Thickness varies sinusoidally with x"""
d_mean = 1.0
amplitude = 0.2
return d_mean + amplitude * np.sin(2 * np.pi * x / Lx)
def thickness_deriv(x):
"""Derivative of thickness: dd/dx"""
amplitude = 0.2
return amplitude * (2 * np.pi / Lx) * np.cos(2 * np.pi * x / Lx)
# Bases - CHANGED to use Z coordinate
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, 1), 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)
# ADDED: Thickness fields
# Create Dedalus fields for thickness and Z coordinate to avoid numpy array issues
d_field = dist.Field(name='d_field', bases=(xbasis, Zbasis))
ddx_field = dist.Field(name='ddx_field', bases=(xbasis, Zbasis))
Z_field = dist.Field(name='Z_field', bases=(xbasis, Zbasis))
# Fill the fields
d_field['g'] = thickness(x)
ddx_field['g'] = thickness_deriv(x)
Z_field['g'] = Z
# Coordinate transformation operators
lift_basis = Zbasis.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
# MODIFIED: Transformed gradient operators
# Original: grad = dx*ex + dz*ez
# Transformed: grad_phys = [dx - (Z/d)*ddx*dZ]*ex + (1/d)*dZ*eZ
grad_u = d3.grad(u) + eZ * lift(tau_u1)
grad_b = d3.grad(b) + eZ * lift(tau_b1)
# Apply transformation to gradients
# For simplicity, we'll construct the transformed operators term by term
dx_op = lambda f: d3.Differentiate(f, coords['x'])
dZ_op = lambda f: d3.Differentiate(f, coords['Z'])
# Transformed buoyancy gradient
grad_b_phys_x = dx_op(b) - Z_field * ddx_field * (1/d_field) * (dZ_op(b) + lift(tau_b1))
grad_b_phys_z = (1/d_field) * (dZ_op(b) + lift(tau_b1))
# Velocity components
ux = u @ ex
uz = u @ eZ
# Transformed velocity gradients
grad_ux_x = dx_op(ux) - Z_field * ddx_field * (1/d_field) * (dZ_op(ux) + lift(tau_u1 @ ex))
grad_ux_z = (1/d_field) * (dZ_op(ux) + lift(tau_u1 @ ex))
grad_uz_x = dx_op(uz) - Z_field * ddx_field * (1/d_field) * (dZ_op(uz) + lift(tau_u1 @ eZ))
grad_uz_z = (1/d_field) * (dZ_op(uz) + lift(tau_u1 @ eZ))
# Divergence and Laplacian
div_u_phys = grad_ux_x + grad_uz_z
laplacian_b = dx_op(grad_b_phys_x) + dZ_op(grad_b_phys_z)
# Problem
problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())
problem.add_equation("div_u_phys + tau_p = 0")
problem.add_equation("dt(b) - kappa*laplacian_b + lift(tau_b2) = - ux*grad_b_phys_x - uz*grad_b_phys_z")
problem.add_equation("dt(u) - nu*(dx_op(grad_ux_x) + dZ_op(grad_ux_z))*ex - nu*(dx_op(grad_uz_x) + dZ_op(grad_uz_z))*eZ + dx_op(p)*ex + (1/d_field)*dZ_op(p)*eZ - b*eZ + lift(tau_u2) = - ux*(grad_ux_x*ex + grad_uz_x*eZ) - uz*(grad_ux_z*ex + grad_uz_z*eZ)")
# Boundary conditions at Z=0 and Z=1
# Need 1D field for BC at bottom
d_1d = dist.Field(name='d_1d', bases=xbasis)
x_1d = dist.local_grid(xbasis)
d_1d['g'] = thickness(x_1d)
problem.add_equation("b(Z=0) = d_1d")
problem.add_equation("u(Z=0) = 0")
problem.add_equation("b(Z=1) = 0")
problem.add_equation("u(Z=1) = 0")
problem.add_equation("integ(p) = 0")
# 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)
b['g'] *= Z * (1 - Z)
# Add linear background that accounts for thickness
b['g'] += thickness(x) * (1 - Z)
# Analysis
snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.25, max_writes=50)
snapshots.add_task(b, name='buoyancy')
snapshots.add_task(d_field, name='thickness')
snapshots.add_task(ux, name='u_x')
snapshots.add_task(uz, name='u_z')
# CFL
CFL = d3.CFL(solver, initial_dt=max_timestep, cadence=10, safety=0.5, threshold=0.05,
max_change=1.5, min_change=0.5, max_dt=max_timestep)
CFL.add_velocity(u)
# Flow properties
flow = d3.GlobalFlowProperty(solver, cadence=10)
flow.add_property(np.sqrt(u@u)/nu, name='Re')
# Main loop
try:
logger.info('Starting main loop')
while solver.proceed:
timestep = CFL.compute_timestep()
solver.step(timestep)
if (solver.iteration - 1) % 10 == 0:
max_Re = flow.max('Re')
logger.info(f"Iteration={solver.iteration}, Time={solver.sim_time:.3e}, dt={timestep:.3e}, max(Re)={max_Re:.3f}")
except Exception as e:
logger.error(f'Exception raised: {e}, triggering end of main loop.')
raise
finally:
logger.info(f"Final iteration={solver.iteration}, sim_time={solver.sim_time:.3e}")
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Keaton Burns
--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/302aca98-34df-475a-b626-125498f73768n%40googlegroups.com.