LHS operator interp_theta is coupled along direction theta.

[Antonin]

Hi guys,

I have this error when I try to run my code, but I don't understand why, I don't have any interpolation operator in the LHS of my equations. Do you know what can cause this error to occur?

Thanks for your help,

Antonin

[Daniel Lecoanet]

Hi Antonin,

Can you give a simple example of where this occurs?

Thanks,

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-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/72d8e238-9b24-4142-87eb-da558028faa0n%40googlegroups.com.

[Antonin]

Hi,

Here is the part of my code where I write the equations:

# Parameters

Lr, Ly = (4., 2*np.pi)

Prandtl = 1.

Rayleigh = 1e6

# Create bases and domain

theta_basis = de.Fourier('theta', 64, interval=(0, Ltheta), dealias=3/2)

r_basis = de.Chebyshev('r', 64, interval=(0, Lr), dealias=3/2)

domain = de.Domain([r_basis, theta_basis], grid_dtype=np.float64)

# 2D Dynamics

problem = de.IVP(domain, variables=['psi','u','v','T0','dT0','w','dw','N','u0','du0'])

problem.parameters['P'] = (Rayleigh * Prandtl)**(-1/2)

problem.parameters['R'] = (Rayleigh / Prandtl)**(-1/2)

problem.parameters['F'] = F = 1

problem.parameters['E'] = E = 1

problem.parameters['L'] = L = 1

problem.parameters['Tm'] = Tm = 1

problem.add_equation("r**2*dt(w) - r*w - r**2*dr(dw) - dtheta(dtheta(w)) - r**2*R*dtheta(T0) = -r**2*u*dr(w) - r*v*dtheta(w)")

problem.add_equation("r**2*dt(u0) + r**2*N + r**2*E**(-0.5)/L*u0 + u0 - r*dr(u0) - r**2*dr(du0) - dtheta(dtheta(u0)) = 0 ")

problem.add_equation("r**2*dt(T0) - P**(-1)*r*dr(T0) - P**(-1)*r**2*dr(dT0) - P**(-1)*dtheta(dtheta(T0)) = -r**2*u*dr(T0) - r*v*dtheta(T0)")

problem.add_equation("du0 - dr(u0) = 0")

problem.add_equation("dT0 - dr(T0) = 0")

problem.add_equation("dw - dr(w) = 0")

problem.add_equation("r*w - r*dr(v) - v + dtheta(u) = 0")

problem.add_equation("v - u0 + dr(psi) = 0")

problem.add_equation("r*N = r*u*dr(v) + v*dtheta(v) + u*v", condition = "(ntheta == 0)")

problem.add_equation("r*u - dtheta(psi) = 0")

problem.add_bc("right(psi) = 0")

problem.add_bc("right(dr(psi)) = 0")

problem.add_bc("right(u0) = 0")

problem.add_bc("left(u0) = 0")

I try to simulate convection in a rotating cylinder with a horizontal temperature gradient.

Antonin

[Daniel Lecoanet]

Hi,

When you build the domain, you need to list the Chebyshev basis last.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/ced57639-81fe-4953-bbed-1cbb5ba45eabn%40googlegroups.com.

[Chiabrando Nicolas]

Hi,

Thank you for your answer, it fixed our problem. However we still find difficulties to launch the program. We have 10 variables, 12 differential equations (10 if we don't count the equations with the condition ntheta == 0/1), and 7 boundary conditions. When we launch the programm, we have this error :" ValueError: Pencil (1,) has 8 constant equations for 0 constant variables plus 7 differential equations / tau terms."

But when we remove 1 boundary condition we have then : "ValueError: Pencil (0,) has 7 constant equations for 0 constant variables plus 8 differential equations / tau terms." and we don't understand why. The problem seems to come from the differential equations with the conditions on ntheta. To be sure, we  have removed those conditions and this time another error message occured : "Factor is exactly singular"". I have seen on internet a similar problem (https://groups.google.com/g/dedalus-users/c/mCYcalFpnoo) but I didn't understand how to fix it. Do you know what can cause those errors to occur ?

Thank you for your help,

Nicolas

Here is our code : 

"""

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_procs` command can

be used to merge distributed analysis sets from parallel runs, and the

`plot_slices.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 -m dedalus merge_procs snapshots

$ mpiexec -n 4 python3 plot_slices.py snapshots/*.h5

This script can restart the simulation from the last save of the original

output to extend the integration. This requires that the output files from

the original simulation are merged, and the last is symlinked or copied to

`restart.h5`.

To run the original example and the restart, you could use:

$ mpiexec -n 4 python3 rayleigh_benard.py

$ mpiexec -n 4 python3 -m dedalus merge_procs snapshots

$ ln -s snapshots/snapshots_s2.h5 restart.h5

$ mpiexec -n 4 python3 rayleigh_benard.py

The simulations should take a few process-minutes to run.

"""

import numpy as np

from mpi4py import MPI

import time

import pathlib

from dedalus import public as de

from dedalus.extras import flow_tools

import logging

logger = logging.getLogger(__name__)

# Parameters

Lr, Ltheta = (4., 2*np.pi)

Prandtl = 1.

Rayleigh = 1e6

# Create bases and domain

theta_basis = de.Fourier('theta', 64, interval=(0.01, Ltheta), dealias=3/2)

r_basis = de.Chebyshev('r', 64, interval=(0.01, Lr), dealias=3/2)

domain = de.Domain([theta_basis, r_basis], grid_dtype=np.float64)

# 2D Boussinesq hydrodynamics

# problem = de.IVP(domain, variables=['psi','u','v','T0','dT0','w','dw','N','u0','du0'])

problem = de.IVP(domain, variables=['psi','u','v','T0','dT0','w','dw','u0','du0', 'N'])

problem.parameters['P'] = 1#(Rayleigh * Prandtl)**(-1/2)

problem.parameters['R'] = 1 #(Rayleigh / Prandtl)**(-1/2)

problem.parameters['F'] = F = 1

problem.parameters['E'] = E = 1

problem.parameters['L'] = L = 1

problem.parameters['Tm'] = Tm = 1

problem.parameters['L'] = Ltheta

problem.parameters['pi'] = np.pi

#Main equations

problem.add_equation("r**2*dt(w) - r*w - r**2*dr(dw) - dtheta(dtheta(w)) - r**2*R*dtheta(T0) = -r**2*u*dr(w) - r*v*dtheta(w)")

problem.add_equation("r**2*dt(u0) + r**2*N + r**2*E**(-0.5)/L*u0 + u0 - r*dr(u0) - r**2*dr(du0) - dtheta(dtheta(u0)) = 0 ")

problem.add_equation("r**2*dt(T0) - P**(-1)*r*dr(T0) - P**(-1)*r**2*dr(dT0) - P**(-1)*dtheta(dtheta(T0)) = -r**2*u*dr(T0) - r*v*dtheta(T0)")

problem.add_equation("r*w - r*dr(v) - v + dtheta(u) = 0")

#Derivatives

problem.add_equation("du0 - dr(u0) = 0")

problem.add_equation("dT0 - dr(T0) = 0")

problem.add_equation("dw - dr(w) = 0")

#psi

problem.add_equation("r*u - dtheta(psi) = 0")

#utheta

problem.add_equation("v - u0 + dr(psi) = 0", condition = "(ntheta == 0)")

problem.add_equation("u0= 0", condition = "(ntheta != 0)")

#N

problem.add_equation("r*N = r*u*dr(v) + v*dtheta(v) + u*v", condition = "(ntheta == 0)")

problem.add_equation("N = 0", condition = "(ntheta != 0)")

problem.add_bc("right(psi) = 0")

problem.add_bc("right(dr(psi)) = 0")

problem.add_bc("left(psi) = 0")

problem.add_bc("left(dr(psi)) = 0")

problem.add_bc("right(u) = 0")

problem.add_bc("left(u) = 0")

problem.add_bc("left(T0) = 0")

problem.add_bc("right(T0) = 0", condition="(ntheta != 1)")

problem.add_bc("right(T0) = 1", condition="(ntheta == 1)")

# Build solver

solver = problem.build_solver(de.timesteppers.RK222)

logger.info('Solver built')

# Initial conditions or restart

if not pathlib.Path('restart.h5').exists():

# Initial conditions

r, theta = domain.all_grids()

T0 = solver.state['T0']

psi = solver.state['psi']

u0 = solver.state['u0']

# 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]

# Linear background + perturbations damped at walls

# zb, zt = z_basis.interval

# pert = 1e-3 * noise * (zt - z) * (z - zb)

# T0['g'] = F * pert

# b.differentiate('z', out=bz)

# Timestepping and output

dt = 0.125

stop_sim_time = 25

fh_mode = 'overwrite'

else:

# Restart

write, last_dt = solver.load_state('restart.h5', -1)

# Timestepping and output

dt = last_dt

stop_sim_time = 50

fh_mode = 'append'

# Integration parameters

solver.stop_sim_time = stop_sim_time

# Analysis

snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.25, max_writes=50, mode=fh_mode)

snapshots.add_system(solver.state)

# # CFL

CFL = flow_tools.CFL(solver, initial_dt=dt, cadence=10, safety=0.5,

max_change=1.5, min_change=0.5, max_dt=0.125, 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) / R", name='Re')

# # Main loop

try:

logger.info('Starting loop')

start_time = time.time()

while solver.proceed:

dt = CFL.compute_dt()

dt = solver.step(dt)

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))