[solver.evaluator.add_file_handler] RunTime error: Factor is exactly singular
louisalexan...@gmail.com
I'm working with Daniel on penetrative convection; you guessed it with use Dedalus, and I have been very happy with it (good job!). Just 2 days ago, however, I noticed something which Daniel called a bug (so maybe you can do something about it for future release).
Basically, the script below breaks at some point (few hundred iterations) during the analysis (saving data). Here is the output:
=================================
Traceback (most recent call last):
File "main2N.py", line 147, in <module>
dt = solver.step(dt, trim=True)
File "/home/lacouston/dedalus/dedalus/core/solvers.py", line 435, in step
self.timestepper.step(self, dt, wall_time)
File "/home/lacouston/dedalus/dedalus/core/timesteppers.py", line 573, in step
p.LHS_LU[i] = linalg.splu(p.LHS.tocsc(), permc_spec=PERMC_SPEC)
File "/home/lacouston/build/lib/python3.5/site-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 251, in splu
ilu=False, options=_options)
RuntimeError: Factor is exactly singular
================================
Here is what Daniel had to say about it:
"It's the trim=True in the timestepper. You're doing an output every 4e-5, and your max_dt=2e-6, so if you're limited by max_dt, you're doing an output every 20 timesteps. If 20*2e-6 is slightly smaller than 4e-5 (due to floating point arithmetic), then the trim ends up setting a tiny timestep. This then causes dedalus to crash. To solve the issue for now, you can make sure max_dt is not an even divisor of your output timescale."
Alternatively, what I have done since was to multiply my time variables by e6, which changes the floating point arithmetic and solved the issue.
Cheers
Louis
******************** SCRIPT ***************************
"""
Dedalus script for 2D Rayleigh-Benard convection.
This script uses a Fourier basis in the x direction with periodic boundary
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.
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
The simulation should take a few process-minutes to run.
"""
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__)
# Problem parameters
Lx, Lz = (4, 2)
Pr=1
Ra=1.0e7
ratio=1
zint=1 # tentative interface
Nu=25 # expected Nusselt number for Ra=?, Pr=1 (c.f. )
Tsz=-Nu # set system close to thermal equilibrium
# Create bases and domain
nx = 512
nz = 128
x_basis = de.Fourier('x',nx, interval=[0, Lx],dealias=3./2.)
z_basis = de.Chebyshev('z',nz, interval=[0, Lz],dealias=3./2.)
domain = de.Domain([x_basis, z_basis], grid_dtype=np.float64)
# def GenFun for variable thermal exp coef
thexp = domain.new_field()
thexp.set_scales(domain.dealias)
thexp.require_grid_space()
def FGen(*args):
np.place(thexp.data,args[0].data<0,-ratio) # if \theta<-1/2: return ratio=\alpha_s/\alpha_c
np.place(thexp.data,args[0].data>=0,1) # if \theta>=-1/2: return 1
return thexp.data*args[0].data
def GGen(*args,domain=domain,F=FGen):
return de.operators.GeneralFunction(domain,layout='g',func=FGen,args=args)
de.operators.parseables['GGen'] = GGen
# 2D Boussinesq hydrodynamics in DIMENSIONLESS SPACE (T is \theta)
problem = de.IVP(domain, variables=['p','T','u','w','Tz','uz'])
problem.meta['p','T','u','w']['z']['dirichlet'] = True
problem.parameters['Pr'] = Pr
problem.parameters['Ra'] = Ra
problem.parameters['ratio'] = ratio
problem.parameters['Tsz'] = Tsz
problem.add_equation("dt(u) - Pr*(dx(dx(u)) + dz(uz)) + dx(p) = -(u*dx(u) + w*uz)")
problem.add_equation("dt(w) - Pr*(dx(dx(w)) - dx(uz)) + dz(p) = -(u*dx(w) + w*dz(w)) + Pr*Ra*GGen(T)")
problem.add_equation("dt(T) - (dx(dx(T)) + dz(Tz)) = -(u*dx(T) + w*Tz)")
problem.add_equation("dx(u) + dz(w) = 0")
problem.add_equation("Tz - dz(T) = 0")
problem.add_equation("uz - dz(u) = 0")
problem.add_bc("left(T) = 1")
problem.add_bc("left(u) = 0")
problem.add_bc("left(w) = 0")
problem.add_bc("right(T) = Tsz")
problem.add_bc("right(u) = 0")
problem.add_bc("right(w) = 0", condition="(nx != 0)")
problem.add_bc("right(p) = 0", condition="(nx == 0)")
# Build solver
solver = problem.build_solver(de.timesteppers.RK222)
logger.info('Solver built')
# Initial conditions
x = domain.grid(0)
z = domain.grid(1)
T = solver.state['T']
Tz = solver.state['Tz']
# 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]
# Piecewise linear background + perturbations damped at walls
zb, zt = z_basis.interval
pert = 1e-3 * noise * (zt - z) * (z - zb)
T['g'] = np.interp(z,[zb, zint, zt],[1, 0, Tsz])+ pert
T.differentiate('z', out=Tz)
# Initial timestep
dt = 4.0e-7
# Integration parameters
solver.stop_sim_time = 0.2 # 1/5 of a diffusive time scale
solver.stop_wall_time = 72 * 60 * 60.
solver.stop_iteration = np.inf
# Analysis
snapshots = solver.evaluator.add_file_handler('snapshots-out', sim_dt=4e-5, max_writes=50)
snapshots.add_system(solver.state, layout='g')
snapshots.add_task("GGen(T)/T", layout='g', name='Thexp')
snapshots.add_task("Ra*GGen(T)", layout='g', name='buoy')
# CFL
CFL = flow_tools.CFL(solver, initial_dt=dt, cadence=10, safety=0.5,
max_change=1.5, min_change=0.5, max_dt=2e-6, threshold=0.05)
CFL.add_velocities(('u', 'w'))
# Flow properties
flow = flow_tools.GlobalFlowProperty(solver, cadence=10)
flow.add_property("sqrt(u*u + w*w) / Pr", 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=True)
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))
******************** END OF SCRIPT ***************************
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 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/4d6ce62e-be06-4319-9646-a23c3b405204%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.