dt, trim=True makes shallow water equations blow up

[louisalexan...@gmail.com]

Hey all,

I was wondering whether we can solve dissipationless wave equations in Dedalus.
And I got scared because all my attempts (changing CFL, max_dt, grid resolution, time step solver) would just blow up (i.e. give nans) even with the simplest equation (e.g. linear 1d shallow water equation).

I think I found the culprit.
If you output stuff and use
dt = solver.step(dt, trim=True)
then it blows up.
If you output stuff and use
dt = solver.step(dt, trim=False)
then it doesn't blow up :)

So trimming of the time step is the culprit.

- If we set trim=False, does that mean that we get the output at the time step closest after the time we wanted the data?
- Do you know if something can be done so we can use Trim=True and not get nans?
- Do you know what's going on??? For more complicated equations, I never had this problem.
- The equations I'm solving (see script below) are solved entirely implicitly (for now). What role does the CFL play then? (I though it's unconditionally stable when you do an implicit solve;; such that it d be constrained by dt_max only)
- Do you expect additional trickeries with the advection term added on the RHS?

ThankS!
Louis

"""
Dedalus script for 2D SWE

"""

import numpy as np
from mpi4py import MPI
import time

from dedalus import public as de
from dedalus.extras import flow_tools

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

import logging
root = logging.root
for h in root.handlers:
h.setLevel("INFO")

logger = logging.getLogger(__name__)
rank = MPI.COMM_WORLD.rank

# Problem parameters
Lx, Ly = (4000, 1000)
nx = 128
ny = 32

# *****USER-PHY***** #
g = 10
f = 1e-3
k = 2*np.pi/1000
H = f**2/(g*k**2)
om = np.sqrt(f**2+g*H*k**2)
# ****************** #

# *****USER-CFL***** #
dt_max = (2*np.pi/om) / 19
CFLfac = 0.4
start_dt = dt_max
# ****************** #

# *****USER-CPU***** #
cpu_time = 0.5 * 60 * 60
end_time = 20000
num_iter = 10000
check_dt = (2*np.pi/om) / 5
# ****************** #

# *****INITIAL***** #
check_num = 0
# ***************** #

# Create bases and domain
x_basis = de.Fourier( 'x', nx, interval=[0, Lx], dealias=3./2.)
y_basis = de.Chebyshev('y', ny, interval=[0, Ly], dealias=3./2.)
domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64)

# 2D Boussinesq hydrodynamics in DIMENSIONLESS SPACE (T is \theta)
problem = de.IVP(domain, variables=['u','v','h'])

problem.meta['v']['y']['dirichlet'] = True

problem.parameters['g'] = g
problem.parameters['f'] = f
problem.parameters['H'] = H

problem.add_equation("dt(u) - g*dx(h) = 0") #+ f*v -0*(u*dx(u)+v*dy(u))
problem.add_equation("dt(v) - g*dy(h) = 0") # - f*u-0*(u*dx(v)+v*dy(v))
problem.add_equation("dt(h) - H*(dx(u)+dy(v)) = 0") #-0*(u*dx(h)+v*dy(h))

problem.add_bc("right(v) = 0")
problem.add_bc("left(v) = 0")

# Build solver
solver = problem.build_solver(de.timesteppers.SBDF2)
logger.info('Solver built')

# Initial conditions
x = domain.grid(0)
y = domain.grid(1)
h = solver.state['h']
h['g'] = 1e-1*np.sin(k*x)

# Initial timestep
dt = start_dt

# Integration parameters
solver.stop_sim_time = end_time
solver.stop_wall_time = cpu_time
solver.stop_iteration = num_iter

# *******Analysis******* #
# Checkpointing
checkpoints = solver.evaluator.add_file_handler('checkpoints', sim_dt=check_dt, max_writes=100)
checkpoints.add_system(solver.state,layout='g')
print(dt_max,check_dt)

# CFL
CFL = flow_tools.CFL(solver, initial_dt=dt, cadence=10, safety=CFLfac,
max_change=1.5, min_change=0.5, max_dt=dt_max, threshold=0.05)
CFL.add_velocities(('u', 'v'))

# Flow properties
flow = flow_tools.GlobalFlowProperty(solver, cadence=10)
flow.add_property("sqrt(u*u + v*v)", name='Re')

# Main loop
try:
logger.info('Starting loop')
start_time = time.time()
while solver.ok:
dt = CFL.compute_dt()
dt = solver.step(dt, trim=False)
if (solver.iteration-1) % 10 == 0:
logger.info('Iteration: %i, Time: %e, dt: %e' %(solver.iteration, solver.sim_time, dt))
logger.info('Max Re = %f' %flow.max('Re'))
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*domain.dist.comm_cart.size))

[Daniel Lecoanet]

Hi Louis,

Normally if the trim=True causes things to blow it, it's because your timestep size evenly divides the output cadence of one of you outputs.  For example, if you output every 1.0 time units, and your timestep size is 0.1, then after 10 timesteps, you will have reached a time of:

1/10+1/10+1/10+1/10+1/10+1/10+1/10+1/10+1/10+1/10

0.9999999999999999

Because trim=True, your next timestep will be of size

1.1102230246251565e-16

which causes a nan.  In practice this is never a problem for complicated problems where the timestep size is set by a CFL.

If you don't have advection terms, then you don't need a CFL.  However, if you want to resolve waves, you should probably set your timestep size based off the wave speed.

Daniel

--
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-users+unsubscribe@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/d6ae632a-7d83-4d23-a8dd-241a3f707d45%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[louisalexan...@gmail.com]

Hi Daniel,

Yes, I know that problem, but I tried such weird output_time to make sure it wouldn't be a multiplier of the time step that I thought that the issue should be different.

+it seems that my time step is set by the CFL (?!? you can run to see it); which I find weird since as we both said everything is solved implicitly...

I did some more checks and if you set trim=True with
dt_max = (2*np.pi/om) / 20
check_dt = (2*np.pi/om) / (7/3)
it doesn't blow up (with check_dt the sim_dt in checkpointing), but
check_dt = (2*np.pi/om) / (2.3333)
does !
Isn't this crazy ???
5.123 and 11/3 at the denominator of check_dt also give blow ups...

I mean it looks like in most cases it will blow up... which makes it hard to come up with a good random number.

Any insights appreciated...

[Daniel Lecoanet]

Hi Louis,

The script as you sent runs fine on my laptop:

2018-05-30 10:20:32,111 __main__ 0/1 INFO :: Iteration: 6831, Time: 1.999678e+04, dt: 1.648291e+00

2018-05-30 10:20:32,111 __main__ 0/1 INFO :: Max Re = 5.372107

2018-05-30 10:20:32,195 solvers 0/1 INFO :: Simulation stop time reached.

2018-05-30 10:20:32,196 __main__ 0/1 INFO :: Iterations: 6834

2018-05-30 10:20:32,196 __main__ 0/1 INFO :: Sim end time: 20001.454829

2018-05-30 10:20:32,196 __main__ 0/1 INFO :: Run time: 85.55 sec

2018-05-30 10:20:32,197 __main__ 0/1 INFO :: Run time: 0.023764 cpu-hr

Is this different from what you're getting?

As for why the timesteps are being set by the CFL limit...  Well that's because you told it to calculate a timestep size based off the velocities u and v, and then you're taking timesteps based off that!  Just remove the dt = CFL.compute_dt() from your main loop.

Daniel

--
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-users+unsubscribe@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/9403fe61-081c-4417-aa05-8e2d7bb409a0%40googlegroups.com.

[louisalexan...@gmail.com]

Hey,

Uhm I made minor changes so I have different number of iterations (see script below), but this one works because trim=False!

If you make it trim=True it should blow up.

"""

a = 0.1
param = [['g','f','k','H','om','a','Lx','Ly'],[g,f,k,H,om,a,Lx,Ly]]
if rank==0: np.save('param',param)
# ****************** #

problem.add_equation("dt(u) - g*dx(h) -f*v = 0") # -0*(u*dx(u)+v*dy(u))
problem.add_equation("dt(v) - g*dy(h) +f*u = 0") #-0*(u*dx(v)+v*dy(v))

problem.add_equation("dt(h) - H*(dx(u)+dy(v)) = 0") #-0*(u*dx(h)+v*dy(h))

problem.add_bc("right(v) = 0")
problem.add_bc("left(v) = 0")

# Build solver
solver = problem.build_solver(de.timesteppers.SBDF2)
logger.info('Solver built')

# Initial conditions
x = domain.grid(0)
y = domain.grid(1)
h = solver.state['h']

h['g'] = a*np.sin(k*x)