Trouble simulating vortex evolution in geostrophic flow

[Ann Feng]

Hi everyone,

I’m trying to simulate vortex motion in a geostrophic flow. The code runs without errors and the initial condition seems to be set correctly, but the flow field doesn’t appear to evolve over time.

I suspect the issue might lie in the following part of the code (see below). Could anyone take a look and help me figure out what might be going wrong? Thanks a lot for any pointers!

Ann

# Fields

psi = dist.Field(name='psi', bases=(x_basis,y_basis))

tau_p = dist.Field(name='tau_p')

# Substitutions

dx = lambda A: de.Differentiate(A, coords['x'])     # Derivatives

dy = lambda A: de.Differentiate(A, coords['y'])

L = lambda A: de.Laplacian(A)                       # Laplacian

J = lambda A, B: dx(A)*dy(B) - dy(A)*dx(B)          # Jacobian

HD = lambda A: -mu*L(L(L(L(A))))                    # Hyperdiffusion:  ν8 ∇^8 ψ

u = de.Differentiate(-psi, coords['y'])             # Velocity

v = de.Differentiate(psi, coords['x'])

# Governing equation

problem = de.IVP([psi, tau_p], namespace=locals())

#problem.add_equation("dt(L(psi) - psi/Ld/Ld) + beta*dx(psi) - HD(psi) + tau_p = -J(psi, L(psi))")

problem.add_equation("dt(L(psi) - psi/Ld/Ld) + beta*dx(psi) + tau_p = 0")

problem.add_equation("integ(psi) = 0")

# Timestepping

timestepper = de.RK443

# Solver

solver =  problem.build_solver(timestepper)

# Integration parameters

solver.stop_sim_time = 300.*86400.      # run for 600 days

solver.stop_wall_time = np.inf

solver.stop_iteration = np.inf

# Initial condition

x0 = 200e3 #0.0e3

y0 = 100e3 #0.0e3

xm, ym = np.meshgrid(x,y, indexing='ij')

tmp_eta = eta_h_dim * np.exp( -((xm-x0)**2+(ym-y0)**2)/(eta_r_dim*eta_r_dim))     # vortex's height profile is Gaussian (amp = 15 cm, L = eta_r_dim)

psi['g'] = tmp_eta*gravity/f0 # Set the initial condition of psi, 'g' for the grid-space values.

[Jeffrey S. Oishi]

Hi Ann,

There's nothing obviously wrong with this section of code (that I can see at least). Could you send a full example so someone can try to run it? Without that it's very hard to see what's going on.

Jeff

--
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/66b4c970-7561-4608-94c2-9795052ceb35n%40googlegroups.com.

[Ann Feng]

(I think my last reply didn’t go through, so I’m posting it again.)

  Thanks for the reply! Here’s the full code:  

jso...@gmail.com 在 2025年9月24日 星期三晚上7:22:48 [UTC+8] 的信中寫道：

[Ann Feng]

I’ve tested the linear case and used vector operations in the governing equation, but the problem remains — the initial field still doesn’t seem to evolve. Any advice would be greatly appreciated!

Ann

import numpy as np

import matplotlib.pyplot as plt

import h5py

import time

from scipy import special

from dedalus import public as d3

from dedalus.extras import flow_tools

import logging

logger = logging.getLogger(__name__)

# Domain size

Nx, Ny = 512, 256

basinscale = 4096.*1.e3

Lx, Ly = basinscale, basinscale*Ny/Nx

# a reasonable scale for ocean vortices

gravity = 9.81

beta    = 2.09e-11*50  # for 24N (latitude)

f0      = 5.93e-5      # for 24N (latitude)

# Gaussian vortex parameters

eta_h_dim = 0.15          # vortex amplitude (i.e SSH = 15 cm)

# Coordinates and distributor

coords = d3.CartesianCoordinates('x','y')

dist   = d3.Distributor(coords, dtype=np.float64)

# Bases (real Fourier since fields are real in physical space)

x_basis = d3.RealFourier(coords['x'], size=Nx, bounds=(-0.75*Lx, 0.25*Lx), dealias=3/2)

y_basis = d3.RealFourier(coords['y'], size=Ny, bounds=(-0.5*Ly,  0.5*Ly ), dealias=3/2)

x, y = dist.local_grids(x_basis, y_basis)

ex, ey = coords.unit_vector_fields(dist)

# Fields

psi = dist.Field(name='psi', bases=(x_basis,y_basis))

tau_p = dist.Field(name='tau_p')

# Substitutions

dx = lambda A: d3.Differentiate(A, coords['x'])

dy = lambda A: d3.Differentiate(A, coords['y'])                    # Laplacian

J = lambda A, B: dx(A)*dy(B) - dy(A)*dx(B)          # Jacobian

u = d3.Differentiate(-psi, coords['y'])             # Velocity

v = d3.Differentiate(psi, coords['x'])

u_vec = u*ex + v*ey

# Governing equation

problem = d3.IVP([psi, tau_p], namespace=locals())

problem.add_equation("dt(lap(psi)) + beta*ex@grad(psi) + tau_p = 0")

problem.add_equation("integ(psi) = 0")

# Timestepping

timestepper = d3.RK443

# Solver

solver =  problem.build_solver(timestepper)

# Integration parameters

solver.stop_sim_time = 300.*86400.      # run for 600 days

solver.stop_wall_time = np.inf

solver.stop_iteration = np.inf

# Initial condition

x0 = 200e3 #0.0e3

y0 = 100e3 #0.0e3

xm, ym = np.meshgrid(x,y, indexing='ij')

tmp_eta = eta_h_dim * np.cos((xm-0.2)/Lx*2*np.pi*5)

psi['g'] = tmp_eta*gravity/f0 # Set the initial condition of psi, 'g' for the grid-space values.

# Analysis

snapshots = solver.evaluator.add_file_handler('./snapshots', sim_dt=0.1*86400., max_writes=50)   # output every 1.5 days

snapshots.add_task(psi, layout='g', name='psi')

snapshots.add_task(u, layout='g', name='u')

snapshots.add_task(v, layout='g', name='v')

snapshots.add_task(d3.Laplacian(psi), layout='g', name='vor')

# CFL

dt = 500

CFL = d3.CFL(solver, initial_dt=dt, cadence=10, safety=0.5, threshold=0.05,

                     max_change=1.5, min_change=0.5, max_dt=dt)

CFL.add_velocity(u_vec)

# Main

try:

    logger.info('Starting loop')

    start_time = time.time()

    while solver.proceed:

        timestep = CFL.compute_timestep()

        solver.step(timestep)

        # Every 10 iterations, print a progress message

        if (solver.iteration-1) % 100 == 0:

            E = 0.5*np.mean((dx(psi).evaluate()['g'])**2 + (dy(psi).evaluate()['g'])**2)

            F = np.fft.rfft2(psi['g'])

            amp = np.real(F)[5,0]

            phase = np.angle(F[5,0])

            logger.info('Iteration: %i, Day: %e, dt: %e, E=%.3e' %(solver.iteration, solver.sim_time/86400., timestep, E))

            logger.info(phase)

except:

    logger.error('Exception raised, triggering end of main loop.')

    raise

finally:

    end_time = time.time()

    logger.info('Iterations: %i' %solver.iteration)

    logger.info('Sim end time: %f' %solver.sim_time)

    logger.info('Run time: %.2f sec' %(end_time-start_time))

    logger.info('Run time: %f cpu-hr' %((end_time-start_time)/60/60*dist.comm_cart.size))

Ann Feng 在 2025年9月24日 星期三晚上7:59:24 [UTC+8] 的信中寫道：

[Jeffrey S. Oishi]

Hi Ann,

OK, I've run your code and I don't see anything obviously wrong with the Dedalus end of things. I suspect your parameters might be such that you shouldn't expect to see any evolution over the timeframe you're running the simulation over.

The uncommented equation has only beta*dx(psi) in the time derivative equation. Using the fact that your initial condition is a Gaussian, its characteristic length is eta_r_dim, so we can get an approximate timescale of eta_r_dim/beta, which is O(10^14) sec. That's a lot more than 150 days!

Jeff

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/99c634cb-cb1f-4eff-a581-852410f98793n%40googlegroups.com.

[Ann Feng]

Hi Jeff,

Thanks so much for your help! This code is actually a modification of an earlier version written in Dedalus2, which normally runs fine (vortex deformation can be observed within about 100 days, both in linear and nonlinear cases). However, I don’t seem to get the same result with Dedalus3.

I’ll leave the old version of the code below in case you’d like to take a look when you have time.

Again, I really appreciate your time and help!

Ann

jso...@gmail.com 在 2025年9月25日 星期四晚上7:35:28 [UTC+8] 的信中寫道：

[Jeffrey S. Oishi]

Hi Ann,

OK, I see, there does appear to be something subtle going on here. I won't have time to look at this much more today, but it does seem that time evolving the Laplacian of psi is behaving differently than I expect. 

Jeff

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/6f9c6f92-5173-441a-b184-92b89a53c1b9n%40googlegroups.com.

[Jeffrey S. Oishi]

Hi Ann,

I'm not sure why your results are differing from your old d2 code, but it does seem like your problem is related to using MKS units. When I rescaled your problem into something like Ld units, it seems to evolve just fine. 

I suggest non-dimensionalizing and trying again. It looks like you are hitting underflow and then your derivative just becomes zero. 

Jeff

[Ann Feng]

Hi Jeff,

Ok, I’ll give it a try and share the results once I’ve updated them.
Thanks again — your feedback has been really helpful!

Ann

jso...@gmail.com 在 2025年9月30日 星期二凌晨2:05:37 [UTC+8] 的信中寫道：